Submitted to ApJ, June 14, 2005 

Preprint typeset using I4TgX style emulateapj v. 6/22/04 



in 
o 
o 

(N 

> 
O 

m 

m 
> 

00 

a^ 
m 
\o 
o 
vn 
o 



o 



A UNIFIED, MERGER-DRIVEN MODEL FOR THE ORIGIN OF STARBURSTS, QUASARS, THE COSMIC X-RAY 
BACKGROUND, SUPERMASSIVE BLACK HOLES AND GALAXY SPHEROIDS 

Philip F. Hopkins', Lars Hernquist', Thomas J. Cox', Tiziana Di Matteo", Brant Robertson', Volker Springel' 

Submitted to ApJ, June 14, 2005 

ABSTRACT 

We develop an evolutionary model for starbursts, quasars, and spheroidal galaxies in which supermassive 
black holes play a dominant role. In this picture, mergers between gas-rich galaxies drive nuclear inflows of gas, 
producing intense starbursts and feeding the growth of supermassive black holes. During this phase, the black 
hole is heavily obscured (a "buried" quasar), but feedback energy from its growth expels the gas, rendering 
the black hole briefly visible as a bright, optical source (a "visible" quasar), and eventually halting accretion 
(a "dead" quasar). The self-regulated growth of the black hole accounts for the observed correlation between 
black hole mass and stellar velocity dispersion in spheroidal galaxies. We show that the quasar lifetime and 
obscuring column density depend on both the instantaneous and peak luminosities of the quasar, and determine 
this dependence using a large set of simulations of galaxy mergers varying the host galaxy properties, orbital 
geometry, and gas physics. 

We use our fits to the lifetime and column density to deconvolve observed quasar luminosity functions and 
obtain the evolution of the formation rate of quasars with a certain peak luminosity, n(Lpeak,z)- In our model, 
quasars spend extended periods of time at luminosities well below their peaks, and so n(Lpeak,z) has a maxi- 
mum, falling off at both brighter and fainter luminosities, corresponding to the "break" in the observed quasar 
luminosity function. We obtain self-consistent fits to hard and soft X-ray and optical quasar luminosity func- 
tions for a model in which «(Lpeak7z) varies with redshift according to pure peak luminosity evolution. From 
this form for n(Lpeak,z), and our simulation results for the luminosity dependence of the quasar lifetime and 
obscuring column, we are able to reproduce many observable quantities, including: the column density dis- 
tribution of both optical and X-ray selected quasar samples, the luminosity function of broad-line quasars in 
X-ray samples and the broad-line (Type I, Type II) fraction as a function of luminosity, the mass function of 
active black holes, the observed distribution of Eddington ratios at both low and high redshift, the present-day 
mass function of relic, inactive supermassive black holes and total black hole mass density, and the spectrum 
of the cosmic X-ray background. In each case, our predictions agree well with observations, matching them 
to higher precision than previous tunable models for quasar lifetimes and obscuration similarly fit to the lumi- 
nosity function. We provide a library of Monte Carlo realizations of our modeling for comparison with a wide 
range of observations, using various selection criteria. 
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1. INTRODUCTION 

The measurement of anisotropies in the cosmic microwave 
background (e.g. Spergel et al. 2003) combined with observa- 
tions of high redshift supernovae (e.g. Riess et al. 1998, 2000; 
Perlmutter et al. 1999) have established a "standard model" 
for the Universe, in which the energy density is dominated by 
an unknown form driving accelerated cosmic expansion, and 
most of the mass is non-baryonic, in a ratio of roughly 5:1 
to ordinary matter On small scales, it is believed that struc- 
ture formed through gravitational instability. In the currently 
favored cold dark matter (CDM) paradigm, objects grow hi- 
erarchically, with smaller ones forming first and then merg- 
ing into successively larger bodies. As baryons fall into dark 
matter potential wells, the gas is shocked and then cools radia- 
tively to form stars and galaxies, in a "bottom-up" progression 
(White & Rees 1978). 

Even with the many successes of this picture, the processes 
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underlying galaxy formation and evolution are poorly under- 
stood. For example, there has yet to be an ab initio calcula- 
tion, starting from an initial state prescribed by the standard 
model, resulting in a population of objects that reproduces 
observed galaxies. However, from the same initial condi- 
tions, computer simulations have yielded a new, successful 
interpretation of the Lyman-alpha forest in which absorption 
in caused by density fluctuations in the intergalactic medium 
(e.g.Cenetal. 1994; Zhang et al. 1995; Hernquist etal. 1996), 
over many orders of magnitude in column density (e.g. Katz et 
al. 1996a), explicitly related to growth of structure in a CDM 
universe (e.g. Croft et al. 1998, 1999, 2002; McDonald et al. 
2000, 2004; Hui et al. 2001; Viel et al. 2003, 2004). This sug- 
gests that the difficulties with understanding galaxy formation 
and evolution lie not in the initial conditions or with the de- 
scription of dark matter, but rather with the physics that has 
been used to model the baryons. 

Observations have revealed regularities in the structure of 
galaxies that point to some of this "missing" physics. Super- 
massive black holes appear to reside at the centers of most 
galaxies (e.g. Kormendy & Richstone 1995; Richstone et al. 
1998; Kormendy & Gebhardt 2001) and the masses of these 
black holes are correlated with either the mass (Magorrian et 
al. 1998; McLure & Dunlop 2002; Marconi & Hunt 2003) or 
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the velocity dispersion (i.e. the Mbh-c relation: Ferrarese & 
Men-itt 2000; Gebhardt et al. 2000; Tremaine et al. 2002) of 
spheroids, demonstrating a direct link between the origin of 
galaxies and supermassive black holes. Simulations which 
follow the self-regulated growth of black holes in galaxy 
mergers (Di Matteo et al. 2005; Springel et al. 2005a) have 
shown that the energy released through this process has a 
global impact on the structure of the merger remnant. If this 
conclusion applies to spheroid formation in general, the sim- 
ulations demonstrate that models for the origin and evolution 
of galaxies must account for black hole growth and feedback 
in a fully self-consistent manner. 

Analytical and semi-analytical modeling ("Silk & ReesI 
[1998; Fabian 1999; Wvithe & Loeb 2002, 2003; 
iBeeehnan & Nath 2005) suggests that, beyond a certain 
threshold, feedback energy from black holes can expel gas 
from the centers of galaxies, shutting down accretion onto 
them and limiting their masses. However, these calculations 
usually ignore the impact of this process on star formation 
and therefore do not explain the link between black hole 
growth and spheroid formation, and furthermore make sim- 
plifying assumptions about the dynamics of such accretion. 
For example, the duration of black hole growth is a free 
parameter, which is fixed either using observational estimates 
or assumed to be similar to e.g. the dynamical time of the 
host galaxy or the e-folding time for Eddington-limited 
black hole growth ts = Mbh/M = 4.5 x 10^/"' (e,/0.1)yr for 
accretion with radiative efficiency er = L/Mc^ ^ 0.1 and 
I = L/Lndd ^ 1 (Salpeter 1964). Moreover, these studies have 
adopted idealized models for quasar light curves, usually 
corresponding to growth at a constant Eddington ratio or 
on-off, "light bulb," scenarios. As we discuss below, less 
restrictive modeling suggests that this phase is actually more 
complex. 

Efforts to mod el quasar accretion and feedback mo re self- 
consistently ('e.g.. lCiotti& Osti'iker 1997. 2001; Granato et aP 
l2004h by treating the hydrodynamical response of gas to 
black hole growth have generally been restricted to idealized 
geometries, such as spherical symmetry, employing simple 
models for star formation and galaxy-scale quasar fueling. 
However, these works have made it possible to estimate duty 
cycles of quasars and shown that the objects left behind have 
characteristics similar to those observed, with quasar feed- 
back being a critical element in reproducing these features 
(e.g. Sazonov et al. 2005; Kawata & Gibson 2005; Cirasuolo 
et al. 2005; for a review, see Ostriker & Ciotti 2005). 

Springel et al. (2005b) have incorporated black hole growth 
and feedback into simulations of galaxy mergers and included 
a multiphase model for star formation and pressurization 
of the interstellar gas by supernovae ( Scringel & Hernaui^d 
l2003h to examine implications of these processes for galaxy 
formation and evolution. Di Matteo et al. (2005) and Springel 
et al. (2005a,b) have shown that gas inflows excited by grav- 
itational torques during a merger both trigger starbursts and 
fuel rapid black hole growth. The growth of the black hole 
is determined by the gas supply and terminates as gas is ex- 
pelled by feedback, halting accretion, leaving a dead quasar 
in an ordinary galaxy. The self-regulated nature of black 
hole growth in mergers explains observed correlations be- 
tween black hole mass and properties of normal galaxies 
(|pi Matteo et al. 2005), as well as the color distribution of el- 
lipticals (Springel et al. 2005a). These results lend support to 
the view that mergers have played an important role in struc- 
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Fig. 1 . — Schematic representation of a "cosmic cycle" for galaxy forma- 
tion and evolution regulated by black hole growth in mergers. 



turing galaxies, as advocated especially by Toomre & Toomre 
(1972) and Toomre (1977). (For reviews, see, e.g., Barnes & 
Hernquist 1992; Barnes 1998; Schweizer 1998.) 

Subsequent analysis by Hopkins et al. (2005a,b,c,d) 
has shown that the merger simulations can account for 
quasar phenomena as a phase of black hole growth. Unlike 
what has been assumed in e.g. semi-analytical studies of 
quasars, the simulations predict complicated evolution for 
quasar lifetimes, fueling rates for black hole accretion, 
obscuration, and quasar light curves. The light curves were 
studied by Hopkins et al. (2005a b), who showed that the 
self-termination process gives observable lifetimes ^ 10^ yr 
for bright optical quasars and predicts a large population 
of obscured sources as a natural stage of quasar evolution, 
as implied by observations (for a review, see Brandt & 
Hasinger 2005). Hopkins et al. (2005b) analyzed simulations 
over a range of galaxy masses and found that the quasar 
light curves and lifetimes are always qualitatively similar, 
with both the intrinsic and observed quasar lifetimes being 
decreasing functions of luminosity, with longer lifetimes at 
all luminosities for higher-mass (higher peak luminosity) 
systems. The dependence of the lifetime on luminosity led 
Hopkins et al. (2005c) to suggest a new interpretation of the 
quasar luminosity function, in which the steep bright-end 
consists of quasars radiating near the Eddington limit and 
is directly related to the distribution of intrinsic peak lumi- 
nosities (or final black hole masses) as has been assumed 
previously (e.g.. Small & Blandford 1992; Haiman & Loeg 
1998; Haiman & Menou 2000; Kauffmann & Haehnelt 200(1 
SomerviUe et al.. .2001: TuUy et al.. .2002: Wvithe & Loeb l 
2003| IVolonteri et alJ l2003t iHaiman. Ouataert. & Bowed 
I03 ICrotorretan]2005l) . but where the shallow, faint-end 
of the luminosity function describes black holes growing 
towards or declining from peak phases of quasar activity, 
with Eddington ratios generally between / ~ 0.01 and 1. The 
"break" in the luminosity function corresponds directly to 
the peak in the distribution of intrinsic quasar properties. As 
argued by Hopkins et al. (2005c d) this new interpretation 
of the luminosity function can self-consistently explain 
various properties of both the quasar and galaxy populations, 
connecting the origin of galaxy spheroids, supermassive 
black holes, and quasars. 

Motivated by these results, and earlier work by many others 
which we summarize below, in this paper we consider a pic- 
ture for galaxy formation and evolution, illustrated schemat- 
ically as a "cosmic cycle" in Figure ^ in which starbursts. 
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quasars, and the simultaneous formation of spheroids and su- 
permassive black holes represent connected phases in the lives 
of galaxies. Mergers are expected to occur regularly in a hi- 
erarchical universe, particularly at high redshifts. Those be- 
tween gas-rich galaxies drive nuclear inflows of gas, trigger- 
ing starbursts and fueling the growth of supermassive black 
holes. During most of this phase, quasar activity is obscured, 
but once a black hole dominates the energetics of the central 
region, feedback expels gas and dust, making the black hole 
visible briefly as a bright quasar Eventually, as the gas is 
further heated and expelled, quasar activity can no longer be 
maintained and the merger remnant relaxes to a normal galaxy 
with a spheroid and a supermassive black hole. In some cases, 
depending on the gas content of the progenitors, the remnant 
may also have a disk (Springel & Hernquist 2005; Robert- 
son et al. 2005a). The remnant will then evolve passively and 
would be available as a seed to repeat the above cycle. As the 
Universe evolves and more gas is consumed, the mergers in- 
volving gas-rich galaxies will shift towards lower masses, ex- 
plaining the decline in the population of the brightest quasars 
from z ^ 2 to the present, and the remnants that are gas-poor 
will redden quickly owing to the termination of star formation 
by black hole feedback (Springel et al. 2005a), so that they re- 
semble elliptical galaxies, surrounded by hot X-ray emitting 
halos (e.g. Cox et al. 2005). 

There is considerable observational support for this sce- 
nario, which has led the development of this picture for the 
co-evolution of galaxies and quasars over recent decades. In- 
frared (IR) luminous galaxies are thought to be powered in 
part by starbursts (e.g. Soifer et al. 1984a,b; Sanders et al. 
1986, 1988a,b; for a review, see e.g. Soifer et al. 1987), 
and the most intense examples locally, ultraluminous infrared 
galaxies (ULIRGs), are invariably associated with mergers 
(e.g. Allen et al. 1985; Joseph & Wright 1985; Armus et al. 
1987; Kleinmann et al. 1988; Melnick & Mirabel 1990; for 
reviews, see Sanders & Mirabel 1996 and Jogee 2004). Ra- 
dio observations show that ULIRGs have large, central con- 
centrations of dense gas (e.g. Scoville et al. 1986; Sargent et 
al. 1987, 1989), providing a fuel supply to feed black hole 
growth. Indeed, some ULIRGs have "warm" IR spectral en- 
ergy distributions (SEDs), suggesting that they harbor buried 
quasars (e.g. Sanders et al. 1988c), an interpretation strength- 
ened by X-ray observations demonstrating the prese nce of 
two n on-thermal point sources in NGC6240 ( Komoss a et al.l 
l2003h . which are thought to be supermassive black holes that 
are heavily obscured at visual wavelengths (e.g. Gerssen et 
al. 2004; Max et al. 2005, Alexander et al. 2005a,b). These 
lines of evidence, together with the overlap between bolomet- 
ric luminosities of ULIRGs and quasars, indicate that quasars 
are the descendents of an infrared luminous phase of galaxy 
evolution caused by mergers (Sanders et al. 1988a), an in- 
terpretation supported by observations of quasar hosts (e.g. 
Stockton 1978; Heckman et al. 1984; Stockton & MacKenty 
1987; Stockton & Ridgway 1991; Hutchings & Neff 1992; 
Bahcall et al. 1994, 1995, 1997; CanaUzo & Stockton 2001). 

However, many of the physical processes that connect the 
phases of evolution in Figure[2are not well understood. Early 
simulations showed that mergers produce objects resembling 
galaxy spheroids (e.g. Barnes 1988, 1992; Hernquist 1992, 
1993a) and that if the progenitors are gas-rich, gravitational 
torques funnel gas to the center of the remnant (e.g. Barnes & 
Hernquist 1991, 1996), producing a starburst (e.g. Mihos & 
Hernquist 1996), but these works did not explore the relation- 
ship of these events to black hole growth and quasar activity. 



While a combination of arguments based on time variability 
and energetics suggests that quasars are produced by the ac- 
cretion of gas onto supermassive black holes in the centers 
of galaxies (e.g. Salpeter 1964; Zel'dovich & Novikov 1964; 
Lynden-Bell 1969), the mechanism that provides the trigger to 
fuel quasars therefore remains uncertain. Furthermore, there 
have been no comprehensive models that describe the tran- 
sition between ULIRGs and quasars that can simultaneously 
account for observed correlations like the Mbh-c relation. 

Here, we study these relationships using numerical simu- 
lations of galaxy mergers that account for the consequences 
of black hole growth. In our simulations, black holes ac- 
crete and grow throughout a merger event, producing com- 
plex, time-varying quasar activity. Quasars reach a peak lumi- 
nosity, Lpeak, during the "blowout" phase of evolution where 
feedback energy from black hole growth begins to drive away 
the gas, eventually slowing accretion. Prior to and following 
this brief period of peak activity, quasars radiate at instan- 
taneous luminosities, L, with L < Lpeak- However, we show 
that even with this complex behavior, the global character- 
istics that determine the observed properties of quasars, i.e. 
lifetimes, light curves, and obscuration, can be expressed as 
functions of L and Lpeak, allowing us to make predictions for 
quasar populations that agree well with observations, support- 
ing the scenario sketched in Figure^ 

In § 121 we discuss our methodology and show how the 
quasar lifetimes and obscuration from our simulations can be 
expressed as functions of the instantaneous and peak lumi- 
nosities of quasars. We also define a set of commonly adopted 
models for the quasar lifetime and obscuration against which 
we compare our predictions throughout. In §|3] we apply our 
models to the quasar luminosity function, using the observed 
luminosity function to determine the distribution of quasar 
peak luminosities, and show that this allows us to simultane- 
ously reproduce the hard X-ray, soft X-ray, and optical quasar 
luminosity functions at all redshifts z < 3, and the distribution 
of column densities in both optical and X-ray samples. In §|3 
we determine the time in our simulations when quasars will 
be observable as broad-line objects, and use this to predict the 
broad-line luminosity function and fraction of broad-line ob- 
jects in quasar samples, as a function of luminosity, as well as 
the mass function of low-redshift, active broad-line quasars. 
In §|5] we estimate the distribution of Eddington ratios in our 
simulations as a function of luminosity, and infer Eddington 
ratios in observed samples at different redshifts. In § |6| we 
use our modeling to predict both the mass distribution and to- 
tal density of present-day relic supermassive black holes, and 
describe their evolution with redshift. In § Q we similarly 
apply this model to predict the integrated cosmic X-ray back- 
ground spectrum, accounting for the observed spectrum from 
^ 1 - lOOkeV. In §|8j we discuss the primary qualitative im- 
plications of our results and propose falsifiable tests of our 
picture. Finally, in § |9j we conclude and suggest directions 
for future work. 

Throughout, we adopt a fl^ = 0.3, fl\ = 0.7, Ho = 
70kms~'Mpc"' (h = 0.7) cosmology. 

2. THE MODEL: METHODOLOGY 

Our model of quasar evolution has several elements, which 
we summarize here and describe in greater detail below. 

• In what follows, a "quasar" is taken to mean the course 
of black hole activity in a single merger event. We use 
the term "quasar lifetime" to refer to the time spent by 
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such a quasar at a given luminosity or fraction of the 
quasar peak luminosity, integrated over all black hole 
activity in a single merger event. This is not meant 
to suggest that this would constitute the entire accre- 
tion history of a black hole - a given black hole may 
have multiple "lifetimes" triggered by different merg- 
ers, with each merger in principle fueling a distinct 
"quasar" with its own lifetime. There is no a priori lu- 
minosity threshold for quasar activity - the time history 
can include various epochs at low luminosities and ac- 
cretion rates. 

We model the galaxy mergers using hydrodynamical 
simulations, varying the orbital parameters of the en- 
counter, the internal properties of the merging galax- 
ies, prescriptions for the gas physics, initial "seed" 
black hole masses of the merging systems, and nu- 
merical resolution of the simulations. The black hole 
accretion rate is determined from the surrounding gas 
(smoothed over the scale of our spatial resolution, 
reaching 20pc in the best cases), i.e. the density and 
sound speed of the gas, and its motion relative to 
the black hole, using Eddington-limited, Bondi-Hoyle- 
Lyttleton accretion theory. The black hole radiates 
with a canonical efficiency = 0.1 corresponding to 
a standard Shakura & Sunvaev (1973) thin disk, and 
we assume that ^ 5% of this radiated luminosity is 
deposited as thermal energy into the surrounding gas, 
weighted by the SPH smoothing kernel (which has 
a profile) over the scale of the spatial resolu- 

tion. This scale is such that we cannot resolve the 
complex accretion flow immediately around the black 
hole, but we adopt this prescription because: (1) it re- 
produces the observed slope and normalization in the 
A^BH-cr relation (Di Matteo et al. 2005), (2) it fol- 
lows from observations, based on estimates of the en- 
ergy contained in highly-absorbed UV portion of the 
quasar SED (e.g., Elvis et al. 1994; Telfer et al. 2002!), 
(3) it follows from theoretical considerations of mo- 
mentum couplin g to dust grains in the dense gas very 
near the quasar (iMurrav et al.ll2005l) and hydrodynam- 
ical simulations of small-scale radiative heating from 
quasar accretion (Ciotti & Ostriker 2001), and (4) even 
if the feedback is initially highly collimated, a driven 
wind or shock in a dense region such as the cen- 
ter of the merging galaxies will rapidly isotropize, so 
long as it is decelerated by gravity and the surround- 
ing medium, allowing the high sound speed within the 
shoc k to equal ize angle-dependent pressure variations 
(e.g.. lKoo & McKee 1990 ). and furthermore initial lo- 
cal distortions will be washed away in favor of triax- 
ial structure determined by the large-scale density gra- 
dients (Bisnovatvi- Kogan & Silich 1991.) . as occurs in 
our simulations. 

For each of our merger simulations, we compute the 
bolometric black hole luminosity and column density 
along ^ 1000 lines of sight to the black hole(s) (evenly 
spaced in solid angle), as a function of time from the 
beginning of the simulation until the system has relaxed 
for ~ 1 Gyr after the merger 

We bin different merger simulations by Lpeak, the peak 
bolometric luminosity of the black hole in the sim- 
ulation, and the conditional distributions of luminos- 



ity, P(L|Lpeak), and column density, P(Nii\L, Lpeak), are 
computed using all simulations that fall into a given bin 
in Lpeak- The final black hole mass (black hole mass 
at the end of the individual merger - subsequent merg- 
ers and quasar episodes could further increase the black 

hole mass) is approximately Mg^ « AfEdd(ipeak) (but 
not exactly, see S I2.4L so we obtain similar results if we 
bin instead by M^^. Our calculation of M^^(Lpeak) al- 
lows us to express our conditional distributions of lumi- 
nosity and column density in terms of either peak lumi- 
nosity or final black hole mass. Critically, we find that 
expressed in terms of Lpeak or Mgj^, there is no system- 
atic dependence in the quasar evolution on the varied 
merger simulation properties - this allows us to calcu- 
late a large number of observables in terms of Lpeak or 

Mgjj without the large systematic uncertainties inher- 
ent in attempting to directly estimate e.g. quasar light 
curves in terms of host galaxy mass, gas fraction, multi- 
phase pressurization of the interstellar medium, orbital 
parameters and merger stage, and other variables. 

• The observed quasar luminosity function is the convo- 
lution of the time a given quasar spends at some ob- 
served luminosity with the rate at which such quasars 
are created. Knowing the distributions P(L|Lpeak) and 
P(A^h|L, Lpeak), we Can calculate the time spent by a 
quasar with some Lpeak at an observed luminosity in 
a given waveband. We use this to fit to observational 
estimates of the bolometric quasar luminosity func- 
tion (/)(L), de-convolving these quantities to determine 
the function «(Lpeak); i-e. the rate at which quasars of 
a given peak luminosity must be created or activated 
(triggered in mergers) in order to reproduce the ob- 
served bolometric luminosity function. 

• Given these inputs, we determine the joint distribu- 
tion in instantaneous luminosity and black hole mass, 
column density distribution, peak luminosity and fi- 
nal black hole mass, as a function of redshift, i.e. 

n(L, Ljy, Mbh, A^h, Lpeak, A/gH I ^t redshifts where 
the observed quasar luminosity function can provide 
the necessary constraint. From this joint distribution, 
we can compute, for example, luminosity functions in 
other wavebands, conditional column density distribu- 
tions, active black hole mass functions and Eddington 
ratio distributions, and relic black hole mass functions 
and cosmic backgrounds. We can compare each of 
these results to those determined using simpler models 
for either the qua sar lifetime or column density distri- 
butions; in § 12.51 we describe a canonical set of such 
models, to which we compare throughout this paper 

2.1. The Simulations 

The simulations were performed with GADGET-2 ("SprmgeQ 
2005), a new version of the paraUe l TreeSPH code GADGET 
(Springel, Yoshida, & White 2001). GADGET-2 is based on 
a fully conservative formulation ( Spr ingel & HernauisL2003) 
of smoothed particle hydrodynamics (SPH), which maintains 
simultaneous energy and entropy conservation when smooth- 
ing lengths evolve adaptively (for a discussion, see e.g., Hern- 
quist 1993b, O'Shea et al. 2005). Our simulations account for 
radiative cooling, heating by a UV background (as in Katz et 
al. 1996b, Dave et al. 1999), and incorporate a sub-resolution 
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model of a multiphase interstellar medium (ISM) to describe 
star fo rmation and supernova feedback ( Scringel & Hernauist 
Feedback from supernovae is captured in this sub- 
resolution model through an effective equation of state for 
star-forming gas, enabling us to stably evolve disks with ar- 
bitrary gas fractions (see, e.g. Springel et al. 2005b; Robert- 
son et al. 2004). In order to investigate the consequences of 
supernova feedback over a range of conditions, we employ 
the scheme of Springel et al. (2005b), introducing a parame- 
ter qEos to interpolate between an isothermal equation of state 
(^Eos = 0) and the full multiphase equation of state (^eos =1) 
described above. 

Supermassive black holes (BHs) are represented by "sink" 
particles that accrete gas at a rate M estimated using an 
Eddington-limited version of Bondi-Hoyle-Lyttleton accre- 
tion theory (Bondi 1952; Bondi & Hoyle 1944; Hoyle & Lyt- 
tleton 1939). The bolometric luminosity of the black hole is 
Lboi = ^rMc^, where = 0.1 is the radiative efficiency. We 
assume that a small fraction (typically « 5%) of Lboi couples 
dynamically to the surrounding gas, and that this feedback is 
injected into the gas as thermal energy, as described above. 

We have performed several hundred simulations of collid- 
ing galaxies, varying the numerical resolution, the orbit of the 
encounter, the masses and structural properties of the merg- 
ing galaxies, initial gas fractions, halo concentrations, and 
the parameters describing star formation and feedback from 
supernovae and black hole growth. This large set of simu- 
lations allows us to investigate merger evolution for a wide 
range of galaxy properties and to identify any systematic de- 
p endence of our modeling. The galaxy models are described 
in lSprrngel et al.i( 2005b) . and we briefly review their proper- 
ties here. 

The progenitor galaxies in our simulations have virial ve- 
locities Kii- = 80, 113, 160, 226, 320, and 500 kms"'. We con- 
sider cases with gas equation of state parameters ^eos = 0.25 
(moderately pressurized, with a mass-weighted temperature 
of star-forming gas ~ lO^^K) and ^eos =10 (the full, "stiff" 
Springel-Hernquist equation of state, with a mass-weighted 
temperature of star-forming gas ^ lO^K), and initial disk gas 
fractions (by mass) of /gas = 0.2, 0.4, 0.8, and 1 .0. Finally, we 
scale these models with redshift, altering the physical sizes of 
the galaxy components and the dark matte r halo concentration 
in accord with cosmological evolution ( Mo. Mao & WTiit^ 
[1998). Details are provided in Robertson et all (l2005bh . and 
here we consider galaxy models scaled appropriately to re- 
semble galaxies of the same Vvii-,/gas,and ^eos at redshifts 
Zgai = 0, 2, 3, and 6. 

For each simulation, we generate two stable, isolated disk 
galaxies, each w ith an extended dark matter halo with a 
lHernauis3 11990) profile, motivated by cosmological simula- 
tions (e.g. Navarro et al. 1996; Busha et al. 2004) and ob- 
servations of halo properties (e.g. Rines et al. 2002, 2002, 
2003, 2004), an exponential disk of gas and stars, and (option- 
ally) a bulge. The galaxies have masses Mvii- = V^^J{lQGHo) 
for Zgai = 0, with the baryonic disk having a mass fraction 
md = 0.041, the bulge (when present) has a mass fraction nib = 
0.0136, and the rest of the mass is in dark matter typically 
with a concentration parameter c = 9.0. The disk scale-length 
is computed based on an assumed spin parameter A = 0.033, 
chosen to be near the mode in the observed A distribution 
((Vitvitska et al. 2002), and the scale-length of the bulge is set 
to 0.2 times the resulting value. In Hopkins et al. (2005a), we 
describe our analysis of simulation A3, one of our set with 



Vvii- = 160kms"', /gas = 1.0, ^Eos = 1-0, and Zgai = 0, a fidu- 
cial choic e with a rotation curve and mass similar to the Milky 
Way, and'Hopki ns et al.l (|2005b'c d) used a set of simulations 
with the same parameters but varying Vyir = 80, 113, 160, 226, 
and 320 kms"', which we refer to below as runs Al, A2, A3, 
A4, and A5, respectively. 

Typically, each galaxy is initially composed of 168000 dark 
matter halo particles, 8000 bulge particles (when present), 
24000 gas and 24000 stellar disk particles, and one BH parti- 
cle. We vary the numerical resolution, with many of our sim- 
ulations using instead twice as many particles in each galaxy, 
and a subset of simulations with up to 128 times as many par- 
ticles. We vary the initial seed mass of the black hole to iden- 
tify any systematic dependence of our results on this choice. 
In most cases, we choose the seed mass either in accord with 
the observed Mbh-c relation or to be sufficiently small that its 
presence will not have an immediate dynamical effect. Given 
the particle numbers employed, the dark matter, gas, and star 
particles are all of roughly equal mass, and central cusps in the 
dark matter and bulge profiles are reasonably well resolved 
(see Fig 2. in Springel et al. 2005b). The galaxies are then set 
to collide from a zero energy orbit, and we vary the inclina- 
tions of the disks and the pericenter separation. 

A representative example of the behavior of the simulations 
is provided in Figure |2] which shows the time sequence of 
a merger involving two bulge-less progenitor galaxies with 
virial velocities of 160 kms"' and initial gas fractions of 20%. 
During the merger, gas is driven to the galaxy centers by 
gravitational tides, fueling nuclear starbursts and black hole 
growth. The quasar activity is short-lived and peaks twice in 
this merger, both during the first encounter and the final co- 
alescence of the galaxies. To illustrate the bright, optically 
observable phase(s) of quasar activity which we identify be- 
low, we have added nuclear point sources in the center at the 
position(s) of the black hole(s) at times T = 1.03, 1.39 and 
1 .48 Gyr, generating a surface density in correspondence to 
the relative luminosities of stars and quasar at these times. 
At other times, the accretion activity is either obscured or the 
black hole accretion rate is negligible. To make the appear- 
ance of the quasar visually more apparent, we have put a small 
part of its luminosity in "rays" around the quasar. These rays 
are artificial and are only a visual guide. 

2.2. Column Densities & Quasar Attenuation 

From the simulation outputs, we determine the obscuration 
of the black hole as a function of time during a merger by cal- 
culating the column density to a distant observer along many 
lines of sight. Typically, we generate ^ 1000 radial lines-of- 
sight (rays), each with its origin at the black hole location and 
with directions uniformly spaced in solid angle dcos (?d0. For 
each ray, we begin at the origin and calculate and record the 
local gas properties using the SPH formalism and move a dis- 
tance along the ray Ar = rjhsmu where rj <\ and /zsmi is the 
local SPH smoothing length. The process is repeated until 
a ray is sufficiently far from the origin (> 100 kpc) that the 
column has converged. We then integrate the gas properties 
along a particular ray to give the line-of-sight column density 
and mean metallicity. We have varied rj and find empirically 
that gas properties along a ray converge rapidly and change 
smoothly for -q = 0.5 and smaller We similarly vary the num- 
ber of rays and find that the distribution of line-of-sight prop- 
erties converges for > 100 rays. 

From the local gas pr operties, we use the multipha se model 
of the ISM described in lSoringel & Hernauis3 (120031) to deter- 
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Fig. 2. — Time sequence from one of our merger simulations (Vyji = 160 kms"', initial gas fraction 20%). Each panel is 80/7^'kpc on a side and shows the 
simulation time in the upper left corner. Brightness of individual pixels gives the logarithm of the projected stellar mass density, while color hue indicates the 
baryonic gas fraction, from 20% (blue) to less than 5% (red). At T = 1.03, 1.39 and 1.48 Gyr, when the black hole could be seen as an optical quasar, nuclear 
point sources are shown, providing a representation of the relative luminosities of stars and the quasar at these times. 



mine the mass fraction in "hot" (diffuse) and "cold" (molecu- 
lar and HI cloud core) phases of dense gas and, assuming pres- 
sure equilibrium, we obtain the local density of the hot and 
cold phases and their corresponding volume filling factors. 
The resulting values are in rough agreement with those of 
|McKee & Ostriker ( 1977). Given a temperature for the warm, 
partially ionized component of the hot-phase ^ 8000 K, deter- 
mined by pressure equilibrium, we further calculate the neu- 
tral fraction of this gas, typically ^ 0.3-0.5. We denote the 
neutral and total column densities as A^hi and A^h, respectively. 
Using only the hot-phase density allows us to place an effec- 
tive lower limit on the column density along a particular line 



of sight, as it assumes a given ray passes only through the 
diffuse ISM, with > 90% of the mass of the dense ISM con- 
centrated in cold-phase "clumps." Given the small volume 
filling factor (< 0.01) and cross section of cold clouds, we ex- 
pect that the majority of sightlines will pass only through the 
"hot-phase" component. 

Using Lboi = e,.Mc^, we model t he intrinsic quasar con- 
tinuum SED following fMarconi et al. (2004), based on opti- 
cal through hard X- ray observations (e.g., E lvis et alJll994t 
'George et a lJ Il998t Vanden Berk et al. 2001*, 'Perolaetall 
2002; Telf eretalJ 12002: Uedaetal. 2003; Vignali et aT| 
120031) . with a reflection component generated by the 
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PEXRAV model JMagdziarz & ZdziarskJl995h . This yields, 
for example, a B-band luminosity logCLs/L©) = 0.80- 
0.067£ + 0.017£2-0.0023/:^ where C = log(Lboi/io)- 12, 
and we take = 4400 A, but as we model the entire intrin- 
sic SED we can determine the bolometric correction in any 
frequency interval. 

We then use a gas-to-dust ratio to determine the extinction 
along a given line of sight at optical frequencies. Observations 
suggest that the majority of reddened quasars have reddening 
curves similar to that of the Small Magellanic Cloud (SMC; 
Hopkins et al. 2004, Ellison et al. 2005), which has a gas- 
to-dust ratio lower than the Mi lky Way by approx imately the 
same factor as its metallicity ( Bouchet et al. 1985). Hence, 
we consider both a gas-to-dust ratio equal to that of the Milky 
Way, (Ab/A^hi)mw = 8.47 x 10"^^ cm^, and a gas-to-dust ra- 
tio scaled by metallicity, Ab/Nhi = (Z/0.02)(Ab/A^hi)mw- In 
both c ases we use the SMC-like reddening curve of iPeil 
The form of the correction for hard X-ray (2-10 
keV) and soft X-ray (0.5-2 keV) luminosities is similar to 
that of the B-band luminosity. We calculate extinction at X- 
ray frequencies (0.03-1 keV) using the photoelectric absorp- 
tion cross sections of iMorris on & McCammon (1983) and 
non-relativistic Compton scattering cross sections, similarly 
scaled by metallicity. In determining the column density for 
photoelectric X-ray absorption, we ignore the inferred ionized 
fraction of the gas, as it is expected that the inner-shell elec- 
trons which dominate the photoelectric absorption edges will 
be unaffected in the temperature ranges of interest. We do not 
perform a full radiative transfer calculation, and therefore do 
not model scattering or re-processing of radiation by dust in 
the infrared. 

For a full comparison of quasar lifetimes and column den- 
sities obtained varying our calculation of A^h, we refer to 
[Hopkins et al. (2005b) (see their Figures 1, 5, & 6), and note 
their conclusion that, after accounting for clumping of most 
mass in the dense ISM in cold-phase structures, the column 
density does not depend sensitively on our assumptions for 
the small-scale physics of the ISM and obscuration - typi- 
cally, the uncertainties in the resulting quasar lifetime as a 
function of luminosity are a factor 2 at low luminosities 
in the B-band, and smaller in e.g. the hard X-ray. Because 
our determination of the quasar luminosity functions is sim- 
ilar using the hard X-ray data alone or the hard X-ray, soft 
X-ray, and optical data simultaneous ly, th e added uncertain- 
ties in our calculation of n(Lpeak) in § l3.2l below owing to the 
uncertainty in our A^h calculation are small compared to the 
uncertainties owing to degeneracies in the fitting procedure 
and uncertain bolometric corrections. 

2.3. The Nu Distribution as a Function of Luminosity 

Next, we consider the distribution of column densities as a 
function of both the instantaneous and peak quasar luminosi- 
ties. For each simulation, we consider A^h values at all times 
with a given bolometric luminosity L (in some logarithmic 
interval in L), and determine the distribution of column den- 
sities at that L weighted by the total time along all sightlines 
with a given A^h- At each L, we approximate the simulated 
distribution and fit it to a lognormal form. 



P(Nh) = 



1 



exp 



-XoiiNulNH) 
24„ 



(1) 



"blowout" phase, and the quasar sweeps away surrounding 
gas and dust to become optically observable. 

We show the resulting median column density Nu at each 
luminosity L in Figure |3l In the upper left panel, simulations 
with Zgai = are shown in black, those with Zgai = 2 in blue, and 
those with Zgai = 3 in yellow. In the upper right, simulations 
with /gas = 0.4 are shown in black, those with /gas = 0.8 in red. 
In the lower left, simulations with ^eos = 0.25 are shown in 
black, those with q^o% = 1 .0 in green. And in the lower right, 
simulations with Kir = 80, 113, 160, 226, 320, and 500 kms"' 
are shown as black asterisks, purple dots, red diamonds, green 
triangles, yellow squares, and red crosses, respectively. Sim- 
ulations with other values for these parameters (not shown 
for clarity, but see e.g. Hopkins et al. [2005d]) show similar 
trends. 

While the increase in typical A^h values with luminosity ap- 
pears to contradict observations suggesting that the obscured 
fraction decreases with luminosity, this is because the rela- 
tionship shown above is dominated by quasars in growing, 
heavily obscured phases. In these stages, the relationship be- 
tween column density and luminosity is a natural consequence 
of the fact that both are fueled by strong gas flows into the 
central regions of the galaxy - more gas inflow means higher 
luminosities, but also higher column densities. During these 
phases, the lognormal fits to column density as a function 
of instantaneous and peak luminosity presented in this sec- 
tion are reasonable approximations, but they break down in 
the brightest, short-lived stages of merger activity when the 
quasar rapidly heats the surrounding gas and drives a power- 
ful wind, lowering the column density, resulting in a bright, 
optically observable quasar Including in greater detail the ef- 
fects of quasar blowout during the final stages of its growth in 
§13 we find that this modeling actually predicts the observed 
decrease in obscured fraction with luminosity. 

The relationship between A^h and L shows no strong sys- 
tematic dependence on any of the simulation parameters con- 
sidered. At most, there is weak sensitivity to ^eos, in the sense 
that the simulations with ^eos = 1.0 have slightly larger col- 
umn densities at a given luminosity than those with q'eos = 
0.25. We derive an analytical model relating both the ob- 
served column density and quasar luminosity to the inflowing 
mass of gas in Hopkins et al. (2005d), by assuming that while 
it is growing, the black hole mass is proportional to the inflow- 
ing gas mass in the galaxy core (which ultimately produces 
the Magorrian et al. [1998] relation between black hole and 
bulge mass), and assuming Bondi accretion, with obscuration 
along a sightline through this (spherically symmetric) gas in- 
flow. Such a model gives the observed correlation between 
A^H and L, and explains the weak dependence of the column 
density-luminosity relation on the ISM gas equation of state. 
The assumptions above give a relationship of the form 



1 /c.x /cL\ 1/3 



(2) 



This provides a good fit for all but the brightest luminosi- 
ties, where quasar feedback becomes important driving the 



where /o ^ 50 is a dimensionless factor depending on the ra- 
diative efficiency, mean molecular weight, density profile, and 
assumed Mbh - o" relation; mu is the mass of hydrogen; Rc 
the radius of the galaxy core (~ lOOpc); and Cs the effective 
sound speed in the central regions of the galaxy. A q'eos =1.0 
equation of state, with a higher effective temperature, results 
in a factor of w 2 larger sound speed in the dens est regions of 
the galaxy than a ^eos = 0.25 equation of state (fSpringel et alJ 
■2005b ). explaining the weak trend seen. In any event, the de- 
pendence is small compared to the intrinsic scatter for either 
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Fig. 3. — The median fitted total (neutral and ionized) column density Nh at each luminosity L in the snapshots from our series of simulations described in 
§|2| We compare changing concentrations and halo properties with redshift Zgai (upper left), gas fractions /gas (upper right), the equation of state parameter (?eos 
(lower left), and virial velocity Vyir (lower right). At lower right, simulations with Vvii = 80, 1 13, 160, 226, 320, and 500 kms"' are shown as black asterisks, purple 
dots, red diamonds, green triangles, yellow squares, and red crosses, respectively. Other than a possible weak sensitivity to ^eos> the column density distribution 
as a function of luminosity shows no systematic dependence on any of the vaiied simulation parameters. 



equation of state in the value of Nh at a given luminosity, and 
further weakens at high luminosity, so it can be neglected. 
What may appear to be a systematic offset in Nh with Vvii- is 
actually just a tendency for larger Vyir systems to be at higher 
luminosities; there is no significant change in the dependence 
of A^H on L. 

We use our large set of simulations to improve our fits (rel- 
ative to those of Hopkins et al. 2005d) to the A^h distribution 
as a function of instantaneous and peak luminosities. Looking 
at individual simulations, there appears to be a "break" in the 
power-law scaling of Nh with L at L ^ 10" L©. We find that 
the best fit to the median column density Nh is then 



, 0-54 

1022.8 £ 
Nh={ )'^""- ,0.43 



ifL< 1O"L0 



(3) 



Either of th^se two relations provides an acceptable fit to 
the plotted Nh distribution if applied to the entire luminosity 
range {y^ /v « 2.8, 3.2 for the first and second relations, re- 
spectively), but their combination provides a significantly bet- 
ter fit (x^/ V K, \ .5), although it is clear from the large scatter 
in Nh values that any such fit is a rough approximation. De- 
spite the complicated form of this equation, it is, in practice, 
similar to omNn ocL" -'^ fit from previous work and A^// (xL}I^ 
analytical scaling over the range of relevant luminosities, but 
is more accurate by a factor ~ 2-3 at low (< \(f Lq) lumi- 
nosities. For comparison, however, we do consider this sim- 



pler form for NniL) as well as our more accurate fit above in 
our subsequent analysis, and find that it makes little difference 
to most observable quasar properties. At the highest luminosi- 
ties, near the peak luminosities of the brightest quasars, the 
scatter about these fitted median Nh values increases, and as 
noted above the impact of the quasar in expelling surrounding 
gas becomes important and column densities vary rapidly. We 
consider this "blowout" phase in more detail in §|3 

We find that any dependence of (the fitted lognormal 
dispersion) on L or Lpeak is not statistically significant, with 
approximately constant (7^^ ~ 0.4 for individual simulations. 
We similarly find no systematic dependence of a^^ on any of 
our varied simulation parameters. However, it is important to 
note that while the dispersion in A^h for an individual simula- 
tion is ~ 0.4, the dispersion in Nh across all simulations 
at a given luminosity is large, ^ 1 dex. Thus, we fit the effec- 
tive <jn„ at a given luminosity for the distribution of quasars 
and find it is ~ 1-2. Although we have slightly revised 
our fits for greater accuracy at low luminosities, we note that 
this relation is shallower than the relation A^h oc L roughly ex- 
pected if Mbh is constant (L cx p oc A^h) or L cx Mbh always, 
and strongly contrasts with unification models which predict 
static obscuration, or evolutionary models in which A'h is in- 
dependent of L up to some threshold (e.g.. iFabian.l999) . 

2.4. Quasar Lifetimes & Sensitivity to Simulation Parameters 

We define the luminosity-dependent quasar lifetime fg = 
fQ(Lmin) as the time a quasar has a luminosity above a certain 
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Fig. 4. — Integrated intiinsic quasar lifetime above a given reference bolo- 
metric luminosity, tQ(L), as a function of luminosity for simulations with host 
galaxies with total mass (top panel) M^-^i = 0.5-2.0 X lO'^ Mq, and simula- 
tions with final black hole masses (bottom panel) M{^^^ = 0.5-2.0 x 10* Mq 
(i.e. similar peak luminosity Lp^ak ~ IO'^Lq). The simulations cover a 
range in equation of state parameter qEOS, initial disk gas fraction /gas, 
galaxy redshift (for scaling of halo properties) Zg-A, and virial velocities 
Vvir = 113 — 160kms"'. The black line in both cases is for a merger in- 
volving Milky Way-like galaxy models, which we refer to as A3, with 
/gas = 1, ?EOS = 1, Zgal = 0, and Vvir = 160kms"'. 
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Fig. 5. — Fits to the quasar lifetime as a function of luminosity from our 
simulations. Upper left shows the intrinsic, bolometric quasar lifetime fg 
of a set of simulations with Lpcak within a factor of 2 of IO'^Lq, in the 
manner of Figure|3 The black histogram shows the geometric mean of these 
hfetimes, and the black histogram in the lower left shows the differential 
hfetime d?/dlogL from this geometric mean. The black thick line in the 
upper left and red line in the lower left show the best-fit to our analytical form, 
d?/dlogL = exp(— L/Lg). Upper right shows the fitted Tq and resulting 
errors in each peak luminosity (final black hole mass) interval, and the best- 
fit power-law to ?g(Lpeak) (red line). Lower right shows the fitted Lg and 
resulting errors in each peak luminosity (final black hole mass) interval, and 
the best-fit proportionality Lq oc Lpcak (red line). 



reference luminosity L,nin; i.e. the total time the quasar shines 
at L > Lmin- For ease of comparison across frequencies, we 
measure the lifetime in terms of the bolometric luminosity, 
L, rather than e.g. the B-band luminosity. Knowing the dis- 
tribution of column densitie s A^h as a function of luminosity 
and system properties (see § I2.3> . we can then analytically or 
numerically calculate the distribution of observed lifetimes at 
any frequency if we know this intrinsic lifetime. Below ~ 1 
Myr, our estimates of fg become uncertain owing to the ef- 
fects of quasar variability and our inability to resolve the local 
small-scale physics of the ISM, but this is significantly shorter 
than even the most rapid timescales ~ 10 Myr of substantial 
quasar evolution. 

As before, we use our diverse sample of simulations to test 
for systematic effects in our parameterization of the quasar 
lifetime. Figure shows the quasar lifetime as a function of 
reference luminosity L^in for both a set of simulations with 
similar total galaxy mass, Mg^i w IO'^Mq, and similar final 

black hole mass (i.e. similar peak quasar luminosity), Mgj^ w 
10^ Mq. In each case, the simulations cover a range in ^eos, 

/gas; Zgal, and Vvii-. 

As Figure |4] demonstrates, at a given Mgai, there is a wide 
range of lifetimes, with a systematic dependence on several 
quantities. For example, for fixed Mgai, a lower ^eos means 
that the gas is less pressurized and more easily collapses to 
high density, resulting in larger Mgjj and longer lifetimes at 
higher luminosities. Similarly, higher /gas provides more fuel 

for black hole growth at fixed Mgai . However, for a given Mg jj, 
the lifetime tg as a function of Lmin is similar across simula- 
tions and shows no systematic dependence on any of the var- 
ied parameters. We find this for all final black hole masses 
in our simulations, in the range 

^BH ^ 1O''-1O'"M0. We 
have further tested this as a function of resolution, comparing 
with alternate realizations of our fiducial A3 simulation with 
up to 128 times as many particles, and find similar results as 
a function of mL,. 
From Figure ^ it is clear that the final black hole mass or 



peak luminosity is a better variable to use in describing the 
lifetime than the host galaxy mass. The lack of any system- 
atic dependence of either the quasar lifetime or A^H(i,i'peak) 
on host galaxy properties implies that our earlier results (Hop- 
kins et al. 2005 a-d) are reliable and can be applied to a wide 
range of host galaxy properties, redshifts, and luminosities, 
although we refine and expand the various fits of these works 
and their applications herein. Furthermore, the large scatter 
in Iq at a given galaxy mass has important implications for 
the quasar correlation function as a function of luminosity, as 
one cannot associate a single quasar luminosity with hosts of 
a given mass (see Lidz et al. 2005). 

Although the truncated power-laws we have previou sly fit- 
ted to tQ using only the A-series simulations (Hopkins et alJ 
2005b) provide acceptable fits to all our runs, we use our new, 
larger set of simulations to improve the accuracy of the fits 
and average over peculiarities of individual simulations, giv- 
ing a more robust prediction of the lifetime as a function of in- 
stantaneous and peak luminosity. For a given peak luminosity 
Lpsik, we consider simulations with an Lpeak within a factor of 
2, and take the geometric mean of their lifetimes tgiL) (we ig- 
nore any points where fg < 1 Myr, as our calculated lifetimes 
are uncertain below this limit). We can then differentiate this 
numerically to obtain df/dlogL (the time spent in a given log- 
arithmic luminosity interval), and fit some functions to both 
curves simultaneously. Figure|5]illustrates this and shows the 
results of our fitting. We find that both the integrated lifetime 
tgiL) and the differential lifetime df/dlogL are well fitted by 
an exponential, 

d?/dlogL = f^exp[-L/L*], (4) 

where both Iq and Lq are functions of M^^ or Lpeak- The best- 
fit such df /dlogL is shown in the figure as a solid line for sim- 
ulations with Lpeak ^ 10"^ Lq, and agrees well with both the 
numerical derivative df/dlogL (lower left, black histogram) 
and the geometric mean tg(L) (upper left, black histogram). 
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This of course implies 

tQ{L) = t*Q e-^^'-odlogL, (5) 

but we are primarily interested in df /dlogL in our subsequent 
analysis. 

Although our fitted lifetime involves an exponential, it is 
in no way similar to the exponential light curve of con- 
stant Eddington-ratio black hole growth or the model in, e.g., 
iHaiman & Loebl ((1998I) . which give df /dlogL = constant^ 

tS < t*Q. 

Our functional form also has the advantage that, although it 
should formally be truncated with df /dlogL = for L > Lpeak, 
the values in this regime fall off so quickly that we can safely 
use the above fit for all large L. Similarly, at L < 10~'*Lpeak, 
df /dlogL falls below the constant fg to which this equa- 
tion asymptotes. Furthermore, in this regime, the fits above 
begin to differ significantly from those obtained by fitting 
e.g. truncated power-laws or Schechter functions. However, 
these luminosities are well below those we generally consider 
and well below the luminosities where the contribution of a 
quasar with some Lpeak is significant to the observed quanti- 
ties we predict. Moreover, this turndown (i.e. the lower value 
predicted by an exponential as opposed to a power-law or 
Schechter function at low luminosities) is at least in part an 
artifact of the finite simulation duration. The values here are 
also significantly more uncertain, as by these low relative ac- 
cretion rates, the system is likely to be accreting in some low- 
efficiency, ADAF state (e.g. Narayan & Yi 1995), which we 
do not implement directly in our simulations. Rather than in- 
troduce additional uncertainties into our modeling when they 
do not affect our predictions, we adopt these exponential fits 
which are accurate at L > 10"^- 10"^ Lpeak- However, for pur- 
poses where the faint-end behavior of the quasar lifetime is 
important, such as predicting the value and evolution of the 
faint-end quasar luminosity function slope with redshift, a 
more detailed examination of the lifetime at low luminosities 
and relaxation of quasars after the "blowout" phase is neces- 
sary, an d we consider these issues separately in Hopkins et al. 

We also note that in iHopkins et al.l (l2005d) we considered 
several extreme limits to our modeling, neglecting all times 
before the final merger and applying an ADAF correction at 
low accretion rates (taken into account a posteriori by rescal- 
ing the radiative efficiency e, with accretion rate, given the 
assumption that such low accretion rates do not have a large 
dynamical effect on the system regardless of radiative effi- 
ciency), and found that this does not change our results - the 
lifetime at low luminosities may be slightly altered but the 
key qualitative point, that the quasar lifetime increases with 
decreasing luminosity, is robust against a wide range of limits 
designed to decrease the lifetime at low luminosities. 

Figure |5l further shows the fitted fg (upper right) and Lq 
(lower right) as a function of peak quasar luminosity for each 
i'peak- We find that Lq, the luminosity above which the life- 
time rapidly decreases, is proportional to Lpeak, 

L*Q = ttiLpeak, (6) 

with a best fit coefficient ai = 0.20 (solid line). The weak 
dependence of fg on Lpeak is well-described by a power-law. 
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Fig. 6. — Predicted quasar lifetime as a function of luminosity compared to 
that obtained in simulations with and without bulges and with different initial 
seed black hole masses. All simulations shown in this plot are initially iden- 
tical to our fiducial A3 (Milky Way-like) case, but with or without an initial 
stellar bulge and with an initial seed black hole mass as labeled. Diamonds 
show the predicted quasar lifetime 1q, a function of the peak luminosity of 
each simulated quasar, determined from the fits shown in Figure|5] Crosses 
show the lifetime determined directly in the simulations. 



The presence or absence of a stellar bulge in the progen- 
itors can have a significant impact on the quasar light curve 
(Springel et al. 2005b), primarily affecting the strength of the 
strong accretion phase associated with initial passage of the 
merging galaxies (e.g. Mihos & Hernquist 1994). Likewise, 
the seed mass of the simulation black holes could have an 
effect, as black holes with smaller initial masses will spend 
more time growing to large sizes, and more massive black 
holes may be able to shut down early phases of accretion in 
mergers in minor "blowout" events. In Figure|6l we show var- 
ious tests to examine the robustness of our fitted quasar life- 
times to these variations. We have re-run our fiducial Milky 
Way-like A3 simulation both with (right panels) and without 
(left panels) initial stellar bulges in the merging galaxies and 
varying the initial black hole seed masses from IO'^-IO^Mq. 
In each case we compare the lifetime fg determined directly 
from the simulations (crosses) to that predicted from our fits 
above (diamonds), based only on the peak luminosity (final 
black hole mass) of the simulated quasar Again, we find that 
varying these simulation parameters can have a significant ef- 
fect on the final black hole mass, but that the quasar lifetime 
as a function of peak luminosity is a robust quantity, indepen- 
dent of initial black hole mass or the presence or absence of a 
bulge in the quasar host. 

We can integrate the total radiative output of our model 
quasars, 

£rad= / L—— dlogL, (8) 

Ji^^ dlogL 
and using our fitted formulae and L^in ^ Lq we find 

E,,i = L*Qt*Q\oge{\-e-''^-^l''Q). (9) 

Knowing Lrad = ^rM^^^c^, we can compare the final black hole 
mass as a function of peak luminosity to what we would ex- 
pect if the peak luminosity were the Eddington luminosity of 
a black hole with mass Mndd, i^Edd = Er-^EddC^/fs, where f^ is 
the Salpeter time for = 0.1. Equating L^ad = ^rM^yf^ with 
the value calculated in Equation|9] and using the definition of 
the Eddington mass at L = Lpeak and our fitted Lq = a^Lpeak, 
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we obtain 



— 22 — !^ = Q;l — loge ! 

MEdd(ipeak) V f5 / 



1.24/7-, 



(10) 



where fj = (Lpeak/IO'^L©)-^" for the power-law fit to tg. 
For our calculations explicitly involving black hole mass, we 
adopt this conversion unless otherwise noted, as we have per- 
formed our primary calculation (i.e. calculated n(Lpeak)) in 
terms of peak luminosity. Moreover, although this agrees well 
with the black hole masses in our simulations as a function of 
peak luminosity (as it must if the fitted quasar lifetimes are 
accurate), this allows us to smoothly interpolate to the high- 
est black hole masses (- a few X 10''- 10"'Mo), which are 
of particular interest in examining the black hole population 
but for which the number of simulations we have with a given 
final black hole mass drops rapidly. 

This gives explicitly the modifications to the black hole 
mass compared to that inferred from the "light bulb" and 
" cons tant Eddington ratio" models which we outline below in 
§ 12.51 in which quasars shine at constant luminosity or follow 

exponential light curves, and for which Mgjj = MEdd(i'peak)/^ 
where I, the (constant) Eddington ratio, is generally adopted. 
The corrections are small, and therefore most of the black hole 
mass is accumulated in the bright, near-peak quasar phase, 
in good agr eement wit h observational estimates (e.g., Solt^ 
[l982; Yu & Tremaine 2002); we discuss this in greater de- 
tail in §|4]and §|6l Furthermore, the increase of fx with de- 
creasing Lpeak implies that lower-mass quasars accumulate a 
larger fraction of their mass in slower, sub-peak accretion af- 
ter the final merger, while high-mass objects acquire essen- 
tially all their mass in the peak quasar phase. This is seen 
directly in our simulations, and is qualitatively in good agree- 
ment with expectations from simulations and semi-analytical 
models in which the Mbh - o" relation is set by black hole 
feedback in a strong quasar phase. Compared to the assump- 
tion that = MEdd(i'peak)^ this formula introduces a small 
but non-trivial correction in the relic supermassive black hole 
mass function implied by the quasar luminosity function and 
"(i-peak) (see §|6}. 

The predictions of our model for the quasar lifetime 
and evolution can be applied to observations which at- 
tempt to constrain the quasar lifetime from individual 
quasars, for example using the prox i mity e f fect in the Lyg 
forest (Baidik, Duncan, & Ostriker 1988"; "Hai man & Ceiil 
|2i002; Jakobsen et al. 2003; Yu & Lu 2005) and multi- 
epoch observations ( Mar tini & Sch neider 2003). However, 
many observations designed to constrain the quasar life- 
time do so not for individual quasars, but using demo- 
graphic or integral arguments based on the population of 
quasars in some luminosity interval (e.g., Soltan 1982; 
"Haehnelt, Nataraian, & Rees 1998; Yu & Tremaine 2002; 
^ & Lu .2004; Porciani. Magliocchetti. & Norberg 2004 



Grazian et alJl2004l) . Our prediction for these observations 
is similar but slightly more complex, as an observed luminos- 
ity function at a given luminosity will consist of sources with 
different peak luminosities Lpeak, but the same instantaneous 
luminosity, L. Furthermore, the lifetime being probed may be 
either the integrated quasar lifetime above some luminosity 
threshold or the differential lifetime at a particular luminos- 
ity. 

For a given determination of the quasar luminosity func- 
tion using our model for quasar lifetimes and some distribu- 
tion of peak luminosities, we can predict the distribution of 



quasar lifetimes as a function of the observed luminosity in- 
terval. Figure shows an example of such a result, us ing t he 
determination of the luminosity function below in § 13.21 at 
redshift z = 0.5. We consider several bolometric luminosities 
spanning the luminosity function from 10''- 1O''*L0, and for 
each, the distribution of sources (peak luminosities), and the 
corresponding distribution of quasar lifetimes. We show both 
the distribution of integrated quasar lifetimes fg (left panel) 
and the distribution of differential quasar lifetimes df/dlogL 
(right panel). The evolution with redshift is weak, with the 
lifetime increasing by ~ 1 .5 -2 at a given luminosity at z = 2. 
There is furthermore an ambiguity of a factor ~ 2, as some 
of the quasars observed at a given luminosity will only be en- 
tering a peak quasar phase, whereas the lifetimes shown are 
integrated over the whole quasar evolution. This prediction is 
quite different from that of the optical quasar phase which we 
describe below in §|3and in Hopkins et al. ( 2005a), as it con- 
siders only the intrinsic bolometric luminosity, but our model- 
ing and the fits provided above for the bolometric lifetime and 
column density distributions should enable the prediction of 
these quantities, considering attenuation, in any waveband. In 
either case, it is clear that the lifetime distribution for lower- 
luminosity quasars is increasingly more strongly peaked and 
centered around longer lifetimes, in good agreement with the 
limited observational evidence from e.g. Adelberger & Stei- 
del (2005). This is a consequence of the fact that in our 
model quasar lifetimes decrease with increasing luminosity. 
The range spanned in the figure corresponds well to the range 
of quas ar lifetimes imp lied by the observations above and oth- 
ers (e . g. lMartini2004i and references therein). 

2.5. Alternative Models of Quasar Evolution 

Our modeling reproduces at least the observed hard X- 
ray quasar luminosity function by construction, since we use 
the observed quasar luminosity functions to dete rmin e the 
birthrate of quasars of a given Lpeak, n(Lpeak), in § 13.21 It is 
therefore useful to consider in detail the differences in our 
subsequent predictions between various models for the quasar 
lifetime and obscuration, in order to determine to what extent 
these predictions are implied by any model that successfully 
reproduces the observed quasar luminosity function, and to 
what extent they are independent of the observed luminosity 
functions and instead depend on the model of quasar evolu- 
tion adopted. To this end, we define two models for the quasar 
lifetime, and two models for the distribution of quasar column 
densities, combinations of which have been commonly used 
in most previous analyses of quasars. 

For the quasar lifetime, we consider the following two 
cases: 



2oo3 
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"Lifiht-Bulb Model" (e.g., 'Small & Blandford' 799^ 



1 



Haiman. Ouataert. & Bow er 2004). The simplest possi- 
ble model for the quasar light curve, the "feast or famine" 
or "light-bulb" model assumes that quasars have only two 
states, "on" and "off." Quasars turn "on", shine at a fixed 
bolometric luminosity L = Lpeak, defined by a "constant" 

Eddington ratio (i.e. Lpeak = ^^bh) constant quasar 
lifetime fg.LB- Models where quasars hve arbitrarily long 
with slowly evolvi n g mean volume emissivity or mean light 
curve (e.g. Small & Blandford 1992VHai man & M enou 20(X)': 
iKauffmann & Haehnelt 2000) are equivalent to the "light 
bulb" scenario, as they still assume that quasars observed at a 
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Fig. 7. — Predicted distribution (fractional number density per logarithmic interval in lifetime) of quasar lifetimes at different bolometric luminosities, for the 
luminosity function determined in §|3]at z = 0.5. Left panel plots the distiibution of integrated lifetimes fg (time spent over the course of each quasar lifetime 
above the given luminosity). Right panel plots the distribution of differential lifetimes d?/dlogL (time spent by each quasar in a logarithmic interval about the 
given luminosity). 



luminosity L radiate at that approximately constant luminos- 
ity over some universal lifetime fg.LB at a particular redshift. 
We adopt / = 0.3 and fe,LB = 10^ yr, as is commonly assumed 
in theoretical work and suggested by observations (given 
this p r ior) (e.g. lYu & Tre maine"2002*. 'Martini"2004', 'SoltanI 
flgsa lYu & Lul l2004r TP orciani, Magliocchetti, & Norberg 
|2p04; Grazian et al. 2004), and similar to the e-folding time 
of a blac k hole with canonical radiative efficiency e, =0.1 
iSalpeterifT964j) or the dynamical time in a typical galactic 
disk or central regions of the merger. These choices control 
only the normalization of n(Lpeak), and therefore do not affect 
most of our predictions. Where the normalization (i.e. value 
of the constant fg or I) is important, we allow it to vary in 
order to produce the best possible fit to the observations. 

"Exponential (Fixed Eddington Ratio) Model." A some- 
what more physical model of the quasar light curve is ob- 
tained by assuming growth at a constant Eddington ra- 
tio, as is commonly adopted in e.g. semi-analytical mod- 
els which attempt to reproduce quasar luminosity functions 
(e.g. "Kauffmann & Haehnelt 2000; Wyithe & Loeb 2003; 
l^lonteri et al. 2003). In this model, a black hole accretes 
at a fixed Eddington ratio / from an initial mass M, to a final 
mass Mf (or equivalently, a final luminosity L/ = / LEddiM/)), 
and then shuts off. This gives exponential mass and luminos- 
ity growth, and the time spent in any logarithmic luminosity 
bin is constant, 

flff/dlog(L) = fs(ln(10)/0 (11) 



for Li < L < Lf. This is true for any exponential light 
curve; i.e. this model includes cases with an exponen- 
tial decline in quasar luminosity), /(f) cx e^'/'*, such as 
that of Haiman & Loeb ( 1998), with only the normalization 
dt /dlog(L) = f, In(lO) changed, and thus any such model will 
give identical results with correspondingly different normal- 
izations. As with the "light-bulb" model, we are free to 
choose the characteristic Eddington ratio and corresponding 
timescale for this lightcurve, and we adopt / = 0.3 (i.e. ~ 
10** yr) in general. Again, however, we allow the normaliza- 
tion to vary freely where it is important, such that these mod- 
els have the best chance to reproduce the observations. For 
our purposes, models in which this timescale is determined 
by e.g. the galaxy dynamical time and thus are somewhat de- 
pendent on host galaxy mass or redshift are nearly identical 
to this scenario. Further, insofar as the dynamical time in- 
creases weakly with increasing host galaxy mass (as, e.g. for 
a spheroid with Mbh oc Mvk ^ aa^/G, where a is the spheroid 

scale length and Mbh « cr'*, such that fdyn ^ a/cr cx cr cx M^^), 
this produces behavior qualitatively opposite to our predic- 
tions (of increasing lifetime with decreasing instantaneous lu- 
minosity), and yields results which are even more discrepant 
from our predictions and the observations than the constant 
(host-galaxy independent) case. 

A wide variety of "light-bulb" or exponential (constant Ed- 
dington ratio) models are possible, allowing for different dis- 
tributions of typical Eddington ratios and/or quasar lifetimes 
(see e.g. Steed & Weinberg 2003 for an extensive comparison 
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of several classes of such models), but for our purposes they 
are essentially identical insofar as they do not capture the es- 
sential qualitative features of our quasar lifetimes, namely that 
the quasar lifetime depends on both instantaneous and peak 
luminosities, and increases with decreasing instantaneous lu- 
minosity. 

We fit both of the simple models above to the observed 
quasar luminosity functions in the same manner described in 
§13 (i-C- in the same manner as we fit our more complicated 
models of quasar evolution), to determine «(Lpeak)LB for the 
"light-bulb" model and «(Lpeak)Edd for the "fixed Eddington 
ratio" model (see Equations 1151 and [T6l respectively). Thus 
all three models of the quasar light curve, the "light-bulb", 
"fixed Eddington ratio", and our luminosity-dependent life- 
times model produce an essentially identical bolometric lu- 
minosity function. 

We also consider two commonly adopted alternative mod- 
els for the column density distribution and quasar obscuration: 

"Standa rd (Luminosity-Independent) Torus" (e.g. 

lAntonuccil \l993 . This is the canonical obscuration 
model, based on obs ervations of local, low-luminosity 
Seyfert galaxies (e.g., 'Risali ti et alJ 11999^. The column 
density distribution is derived from the torus geometry, where 
we assume the torus inner radius lies at a distance Rj from 
the black hole, with a height Hj, and a density distribution 
p{0) (X exp(-7| cos^l), where 9 is the polar angle and the 
torus lies in the 9 = plane. This results in a column density 
as a function of viewing angle of 

Niii9) =Nh,o exp(-7 1 cos 6* | ) cos(90 - 9) 

.^(^)^.ec=(90-.,((g) = -l) ,12, 

iTreister et alJl2004l) . Here, A^h.o is the column density along 
a line of sight through the torus in the equatorial plane and 7 
parameterizes the exponential decay of density with viewing 
angle. This is a phenomenological model, and as a result the 
parameters are essentially all free. We adopt typical values, an 
equatorial column density A^h.o = lO^'^cm"^, radius-to-height 
ratio /^x/^fr = 11. and density profile 7 = 4. This combination 
of parameters follows Treister et al. (2004), and is designed to 
fit the observed X-ray column density distribution and give a 
ratio of obscured to unobscured quasa rs ~ 3, simila r to the 
mean locally observed value (e.g. Risal iti et alJl999l) . 

"Receding (Luminosity-Dependent) Torus" (e.g. 
iLawrenc^ 1199 Ih . Many observations suggest that the 
fraction of obscure d objects depends on luminosity 
Jlteffen et al. 2003: lUedaetalJ 2003; Hasineer 200l 
lyrimes. Rawlinas. & Willott 2004; S azonov & RevnivtsevI 
12004: Barger et al. 2005; Simpson 20051) . Therefore, some 
theoretical works have adopted a "receding torus" model, 
in which the torus radius Rj (i.e. distance from the quasar) 
is allowed to vary with luminosity, but the height and other 
parameters remain constant. The torus radius is assumed to 
increase with luminosity, enlarging the opening angle and 
thus the fraction of unobscured quasars. In this case, the 
column densities are identical to those shown above, but now 
Rj/Ht = (L/L())" ^ where Lq « 10" Lq is the luminosity at 
which the ratio of obscured to unobscured quasars is w 3 : 1 
and the power-law slope is chosen to fit the dependence of 
obscured fraction on luminosity. 



Both of these column density distributions represent phe- 
nomenological models with several free parameters, explic- 
itly chosen to reproduce the observed differences in quasar 
luminosity functions and column density distributions. De- 
spite this, it is not clear that these functional forms represent 
the best possible fit to the observations they are designed to 
reproduce. Furthermore, comparison of our results in which 
column density distributions depend on luminosity and peak 
luminosity elucidates the importance of proper modeling of 
the dependence of column density on quasar evolution. 

3. THE QUASAR LUMINOSITY FUNCTION 

3.1. The Ejfect of Luminosity-Dependent Quasar Lifetimes 

Given quasar lifetimes as functions of both instantaneous 
and peak luminosities, the observed quasar luminosity func- 
tion (in the absence of selection effects) is a convolution of 
the lifetime with the intrinsic distribution of sources with a 
given Lpeak- If sources of a given L are created at a rate h(L, t) 
(per unit comoving volume) at cosmological time f/f ~ 1 / H(z) 
and live for some lifetime /S.tQ{L), the total comoving number 
density observed will be 

PtH+^tQ{L) 

An= ri(L,t)dt, (13) 

which, for a cosmologically evolving ri(L, f ), can be expanded 
about ri(L,tH), yielding An = liiL^tn) AtqiL) to first order in 
AtQ(L)/tH- Considering a complete distribution of sources 
with some Lpeak, we similarly obtain the luminosity function 

'^^'^^=d^^'^^" / ^^^i^^«(Wdl0g(ipeak)- (14) 

Throughout, we will denote the differential luminosity func- 
tion, i.e. the comoving number density of quasars in some log- 
arithmic luminosity interval, as (/) = d$/dlogL. Here, n(Lpeak) 
is the comoving number density of sources created per unit 
cosmological time per logarithmic interval in Lpeak, at some 
redshift, and df/dlogL is the differential quasar lifetime, i.e. 
the total time that a quasar with a given Lpeak spends in a loga- 
rithmic interval in bolometric luminosity L. This formulation 
implicitly accounts for the "duty cycle" (the fraction of active 
quasars at a given time), which is proportional to the lifetime 
at a given luminosity. Corrections to this formula owing to fi- 
nite lifetimes are of order (df/dlogL)/f//, which for the lumi- 
nosities and redshifts considered here (except for Figure HlT i. 
are never larger than ^1/5 and are generally <C 1, which 
is significantly smaller than the uncertainty in the luminosity 
function itself. 

We next consider the implications of our luminosity- 
dependent quasar lifetimes for the relation between the ob- 
served luminosity function and the distribution of peak lumi- 
nosities (i.e. intrinsic properties of quasar systems). In tradi- 
tional models of quasar lifetimes and light curves, this relation 
is trivial. For example, models in which quasars "turn on" at 
fixed luminosity f or so me fixed lifetime (i.e. the "light-bulb" 
model defined in § 12. 5> imply 

n(Lpeak)LB OC 0(L = Lpeak), (15) 

and models in which quasar light curves are a pure exponen- 
tial growth or decay with some cutoff(s) (e.g., exponential or 
fixed Eddington-ratio models) imply 

dcj) 

n(Lpeak)Edd OC — — . (16) 

dlOg(L) [_^r^ 
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Fig. 8. — We reproduce (thin histogram) the luminosity function of 
lUeda et al. ( 2003) at redshift 7, = 0.5 (thin curve) using the binned differ- 
ential quasar lifetime d</dlogL directly from our simulations and a fitted 
distribution of peak luminosities n(Lp^-^) (thick histogram). For each bin in 
log(Lpea]c), we average the binned differential lifetime of a set of simulations 
with peak luminosity in the bin. This clearly demonstrates our key qualitative 
result, that the faint end of the luminosity function is reproduced by quasars 
with peak luminosity around the break luminosity but observed primarily in 
sub-Eddington states (luminosities L <^ Lpcak). is not an artifact of our fitting 
formulae or extrapolation to extreme luminosities. 



These both have essentially identical shape to the observed lu- 
minosity function, qualitatively different from our model pre- 
diction that n(Lpeak) should turn over at luminosities approx- 
imately below the break in the observed luminosity function 
(see, e.g. Fig. 1 of Hopkins et al. 2005e). The luminosity- 
dependent quasar lifetimes determined from our simulations 
imply a new interpretation of the luminosity function, with 
n(Lpeak) tracing the bright end of the luminosity function sim- 
ilar to traditional models, but then peaking and turning over 
below Lpeak ~ i-break, the break luminosity in standard dou- 
ble power-law luminosity functions. In our deconvolution of 
the luminosity function, the faint end corresponds primarily 
to sources in sub-Eddington phases transitioning into or out 
of the phase(s) of peak quasar activity. There is also some 
contribution to the faint-end lifetime from quasars accreting 
efficiently (i.e. growing exponentially at high Eddington ratio) 
early in their activity and on their way to becoming brighter 
sources, but this becomes an increasingly small fraction of 
the lifetime at lower luminosities. For example, in Figure 
7 of Hopkins et al. (2005b), direct calculation of the quasar 
lifetime shows that sub-Eddington phases begin to dominate 
the lifetime for L < 0.1 Lpeak, with > 90% of the Hfetime at 
L ^ 10"^ Lpeak corresponding to sub-Eddington growth. By 
definition, a "fixed Eddington ratio" or "light bulb" model is 
dominated at all luminosities by a fixed, usually large, Ed- 
dington ratio. Even models which assume an exponential de- 
cline in the quasar luminosity from some peak, although they 
clearly must spend a significant amount of time at low Ed- 
dington ratios, have an identical n(Lpeak) = n(Lpeak)Edd (mod- 
ulo an arbitrary normalization), and predict far less time at 
most observable (> 10"'* Lpeak) low luminosities and accretion 
rates (because the accretion rates fall off so rapidly); i.e. the 
population at any observed luminosity is still dominated by 
objects near their peak. 

From our new, large set of simulations, we test this model 
of the relationship between the distribution of peak quasar lu- 
minosities and observed luminosity functions, namely our as- 
sertion that n'(Lpeak) should peak around the observed break 
in the luminosity function, and turn over below this peak. 



with the observed luminosity function faint-end slope dom- 
inated by sources with peak luminosities near the break in 
sub-Eddington (sub-peak luminosity) states. In particular, we 
wish to ensure that this behavior for n(Lpeak) is real, and not 
some artifact of our fitting functions for the quasar lifetime. 

Figure|8lshows the best fit n(Lneak) distribution (solid thick 
histogram) fitted to the lUeda et al. (2003) hard X-ray quasar 
luminosity function (solid curve) at redshift z = 0.5, as well 
as the resulting best-fit luminosity function (solid thin his- 
togram). For ease of comparison with other quasar lumi- 
nosities, we rescale the luminosity function to the bolomet- 
ric luminosity using the corrections of Marconi et al. (2004). 
We determine n'(Lpeak) by logarithmically binning the range of 
Lpeak, and considering for each bin all simulations with Lpeak 
in the given range. For each bin, then, we take the average 
binned time the simulations spend in each luminosity inter- 
val, and take that to be the quasar lifetime df /dlogL. We then 
fit to the observed luminosity function of lUeda etaDlHool, 
fitting 



»-(^peak,)^ AlogL / ^''^ 



and allowing n,(Lpeak. ;) to be a free coefficient for each binned 
Lpeak = Lpeak. i- Dcspitc our large number of simulations, the 
numerical binning process makes this result noisy, especially 
at the extreme ends of the luminosity function. However, the 
relevant result is clear - the qualitative behavior of n(Lpeak) 
described above is unchanged. For further discussion of the 
qualitative differences between the «(Lpeak) distribution from 
different quasar models, and the robust nature of our inter- 
pretation even under restrictive assumptions (e.g. ignoring the 
early phases of merger activity or applying various models for 
radiati ve efficiency as a function of accretion rate), we refer 
to .HoDkins et alJ j2005d) . 

3.2. The Luminosity Function at Different Frequencies and 
Redshifts 

Given a distribution of peak luminosities n(Lpeak), we can 
use our model of quasar lifetimes and the column density dis- 
tribution as a function of instantaneous and peak luminosi- 
ties to predict the luminosity function at any frequency. From 
a distribution of A^h values and some a priori known mini- 
mum observed luminosity L™", the fraction /obs of quasars 
with a peak luminosity Lpeak and instantaneous bolometric lu- 
minosity L which lie above the luminosity threshold is given 
by the fraction of A^h values below a critical A'JJ'^'', where 
L™n=/^Lexp(-cr^A^™''). Here, /^(L) = L^ /Lis abolometric 
correction and ai, is the cross-section at frequency v. Thus, 



±lnf^ 

CTu V L 



fJL)L\ 

min J ' 



and for the lognormal distribution above. 



1 r 



l+erf| 



-log(jVg^VjV„) 



(18) 



(19) 



This results in a luminosity function (in terms of the bolomet- 
ric luminosity) 



./obs (^1^5 ^peak ; ) 



dlog(L) 



n(Lpeak)'^log(Lpeak), (20) 
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where 0(1/. L,L™") is the number density of sources with 
bolometric luminosity L per logarithmic interval in L, with 
an observed luminosity at frequency v above L™". 

Based on the direct fit for n(Lpeak) in Figure |8] we wish 
to consider a functional form for «(Lpeak) with a well-defined 
peak and falloff in either direction in log(Lpeak)- Therefore, 
we take n(Lpeak) to be a lognormal distribution, with 



n(ipeak) = n* 



1 



'111 



exp 



'log(Lpeak/i*) 



(21) 



Here, n» is the total number of quasars being created or acti- 
vated per unit comoving volume per unit time; is the center 
of the lognormal, the characteristic peak luminosity of quasars 
being born (i.e. the peak luminosity at which «(Lpeak) itself 
peaks), which is directly related to the break luminosity in the 
observed luminosity function; and ci* is the width of the log- 
normal in n'(Lpeak)i and determines the slope of the bright end 
of the luminosity function. Since our model predicts that the 
bright end of the luminosity function is made up primarily of 
sources at high Eddington ratio near their peak luminosity, i.e. 
essentially identical to "light-bulb" or "fixed Eddington ratio" 
models, the bright-end slope is a fitted quantity, determined by 
whatever physical processes regulate the bright-end slope of 
the active black hole mass function (possibly feedback from 
outflows or threshold cooling processes, e.g. Wyithe & Loeb 
2003; Scannapieco & Oh 2004; Dekel & Birnboim 2004), un- 
like the faint-end slope which is a consequence of the quasar 
lifetime itself, and is only weakly dependent on the underly- 
ing faint-end active black hole mass or n(Lpeak) distribution. 

We note that although this choice of fitting function has ap- 
propriate general qualities, it is ultimately somewhat arbitrary, 
and we choose it primarily for its simplicity and its capacity to 
match the data with a minimum of free parameters. We could 
instead, for example, have chosen a double power-law form 
with n(Lpeak) = n*/[(ipeak/i*)'" +(ipeak/i*)'''] and 71 < 72, 
but given that the entire faint end of the luminosity function 
is dominated by objects with Lpeak ~ i*, the observed lumi- 
nosity function has essentially no power to constrain the faint 
end slope 71, other than setting an upper limit 71 < 0. The 
"true" n(Lpeak) will, of course, be a complicated function of 
both halo merger rates at a given redshift and the distribution 
of host galaxy properties including, but not necessarily lim- 
ited to, masses, concentrations, and gas fractions. 

Having chosen a form for «(Lpeak), we can then fit to 
an observed luminosity function to determine {nf.,Lf.,af,). 
We take advantage of the capability of our model to pre- 
dict the luminosity function at multiple frequencies, and con- 
sider both fits to just the l.'eda et al. (2003) hard X- ray (2- 
10 keV ) luminosity f u nction , (pHX, and fits to the Ue da et al.l 
^&), 'M ivaii et alJ HmS soft X-ray (0.5-2 keV; (j^sx), 
and Croom et all ^QA) optical B-band (4400 A; 0b) lumi- 
nosity functions simultaneously. These observations agree 
with other, more recent dete rminations of ^hx. <i> sx, <t>ij 
(e.g.lBargeret alJl2005t iHasin ger. Mivaii, & Schmidt 120051 
lRichards'etani2005[ respectivelv^ at most luminosities, and 
therefore we do not expect revisions to the observed lumi- 
nosity functions to dramatically change our results. In order 
to avoid numerical artifacts from fitting to extrapolated, low- 
luminosity slopes in the analytical forms of these luminosity 
functions, we directly fit to the binned luminosity function 
data. Thus, we fit each luminosity function in all redshift in- 
tervals for which we have binned data. 

We find good fits (x^/v = 68.8/104 « 0.66) to afl lumi- 
nosity functions at all redshifts with a pure peak-luminosity 



evolution (PPLE) model, for which 

L, = L^l exp(kiT), ri^ = constant, cr, = constant, (22) 

where r is the fractional lookback time (r = Hq J^dt) and 
Icl is a dimensionless constant fitted with L^,,ri^,, cr*. It is 
important to distinguish this from "standard" pure luminos- 
ity evolution (PLE) models (e.g., Bovle et al. 1988), as with 
n(Lpeak) > and L» = L*(z) always, the density of sources, 
especially as a function of observed luminosity at some fre- 
quency, evolves in a non-trivial manner. 

We do not find significant improvement in the fits if we ad- 
ditionally allow n» or ci* to evolve with redshift (A^^ ~ 1 - 2, 
depending on the adopted form for the evolution), and there- 
fore consider only the simplest parameterization above (Equa- 
tion |23. We also find acceptable fits for a pure density evo- 
lution model, with L, = constant and = «'° expik^T) (both 
keeping ct, fixed and allowing it to evolve as well). However, 
the fits are somewhat poorer ix^/i' ~ 1), and the resulting 
parameters over-produce the present-day density of low-mass 
supermassive black holes and the intensity of the X-ray back- 
ground by an order of magnitude, so we do not consider them 
further In either case, there is a considerable degeneracy be- 
tween the parameters cr, and L, , where a decrease in can be 
compensated by a corresponding increase in cr, . This degen- 
eracy is present because, as indicated above, the observed lu- 
minosity function only weakly constrains the faint-end slope 

of «'(Lpeak)- 

The observations shown are insufficient at high redshift to 
strongly resolve the "turnover" in the total comoving quasar 
density at z ~ 2 - 3, and thus we acknowledge that there must 
be corrections to this fitted evolution at higher redshift, which 
we address below. However, as we primarily consider low 
redshifts, z < 3, and show that the supermassive black hole 
population and X-ray background are dominated by quasars 
at redshifts for which our «(Lpeak) distribution is well deter- 
mined, this is not a significant source of error in most of our 
calculations even if we extrapolate our evolution to 3. 

Figure |9l shows the resulting best-fit PPLE luminosity 
functions from the best-fit «(Lpeak) distribution, for red- 
shifts z = 0-3. This has the best-fit (y^ jv = 0.67) values 
(logU,kL, logri,, CT*) = (9.94, 5.61,-6.29, 0.91) with corre- 
sponding errors (0.29,0.28,0.13,0.09). Here, is in so- 
lar luminosities and ri^, in comoving Mpc"^ Myr"' . Fit- 
ting to the hard X-ray data alone gives a similar fit, 
with the slightly different values (log L*, log cr*) = 
(9.54,4.90,-5.86, 1.03)±(0.66, 0.43, 0.37, 0.13), xVj^ = 0.7 
(note the degeneracy between and cr* in the two fits). Our 
best-fit value of fc^ = 5.6 compares favorably to the value ~ 6 
found by e.g. Boyle et al. (2000) and Croom et al. (2004) for 
the evolution of the break luminosity in the observed luminos- 
ity function, demonstrating that the break luminosity traces 
the peak in the n(Lpeak) distribution at all redshifts. These fits 
and the errors were obtained by least-squares minimization 
over all data points (comparing each to the predicted curve at 
its redshift and luminosity), assuming the functional form we 
have adopted for n(Lpeak)- 

The agreement we obtain at all redshifts, in each of the hard 
X-ray (black solid line), soft X-ray (red dashed line), and B- 
band (dark blue dotted line) is good. This is not at all guaran- 
teed by our procedure, as the fit is highly over-constrained, be- 
cause we fit three luminosity functions each at five redshifts to 
only four free parameters. Of course, the choice of the func- 
tional form for n(Lpeak) ensures that we should be able to re- 
produce at least one luminosity function and its evolution (e.g. 



16 



Hopkins et al. 




Log( Bolometric Luminosity / L®) 



Fig. 9. — Best-fit luminosity function from the pure peak-luminosity evolution n(Lpeak) distribution, for redshifts z = — 3. From our fitted lognormal n(Lpi;ak) 
distribution, we simultaneously reproduce the luminosity function in the hard X-ray (2-10 keV; solid black line), soft X-ray (0.5-2 keV; dashed red line), and 
optical B-band (4400 A; dotted blue fine) at all redshifts. Moreover, we reproduce the distribution of broad-line quasars in hard X-ray selected samples (cyan 
dot-da shed line), as described in §|4] All quantities have been rescaled to bolometric luminosities for ease of comparison, using the connections of Marconi et a^ 
(2004"), with the plotted error bars representing both quoted measurement errors and the estim ated errors in the bolometric corrections. The observ ations are 
from Mivaii et al. ( 200 1 ) (soft X-ray; red squares), Ueda et al. 1 2003) (hard X-ray; black circles'). FCroom et^Il |2004) (B-band, blue diamonds), and, Barger et alj 
-ray selected broad-line quasars; cyan crosses). 



the hard X-ray luminosity function, which is least affected by 
attenuation), but our modeling of the column density distribu- 
tions in mergers allows us to simultaneously reproduce the lu- 
minosity functions in different wavebands without imposing 
assumptions about obscured fractions or sources of attenua- 
tion. Expressed as bolometric luminosity functions, <j)B, <j>sx, 
and <j)HX would be identical in the absence of obscuration, 
similar to the predicted (j)Hx as obscuration is minimal in the 
hard X-ray. 

For redshifts z < 1 , we reproduce in our Figure^] Fig. 2 
of iHopkins et al.l ( l2005dh . w hich shows in detail the agree- 
ment between har d X-ray JUeda et alJ 1200 .3). soft X-ray 
jMivaii et alJl200 tf). and optical (" Bovle et al.ll2000 ) luminos- 
ity functions resulting from the time and luminosity depen- 
dent column density distributions derived from the simula- 
tions. The differential extinction predicted for different fre- 
quencies (and magnitude limits) of observed samples based 
on the column density distributions in our simulations ac- 
counts for the different shape of the luminosity function in 
each band, and the evolution of the luminosity function with 
redshift is driven by a changing L,, the peak of the «(Lpeak) 
distribution (Equation I22> . We emphasize that in our anal- 
ysis, the key quantity constrained by observations is the fit- 
ted «(Lpeak) distribution with redshift. All other quantities 
and distributions are derived from the basic input physics of 
our simulations, with no further assumptions or adjustable 
factors in our modeling beyond the prescription for Bondi 
(Eddington-limited) accretion and ~ 5% energy deposition in 
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Fig. 10. — Hard X-ray (thick), soft X-ray (thin), and B-band (dot-dash) LFs 
determined from our model of quasar lifetimes and column densities, based 
on a distribution of intrinsic source properties fitted to the observed hard X- 
ray LF and the limiting magnitudes of observed samples, at the different red- 
shifts shown. All quantities are rescaled to bolometric luminosities with the 
bolometric connections of Marconi et al. |2004). Symbols show the ob served 
LFs for hard X-rays lUeda et al. 2003, diamonds), soft X-rays lIMivaii et alj 
2000, triangles), and B-band iBovle et al. 20W, crosses). Reproduced from 
iHopkins et aht2005di) . 



the ISM, which are themselves constrained by observations 
and theory as discussed in §|5]and in Di Matteo et al. (2005). 

We can, of course, fit the previously defined simpler model 
of quasar lifetimes, either a "light-bulb" or exponential light 
curve/fixed Eddington ratio model, and obtain an identi- 
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cal hard X-ray luminosity function. We determine these 
fits (see also Equation [21 & I16> and use them through- 
out when we com pare the predictions of such models (de- 
scribed in § 12. 5> to those of our simulated quasar life- 
times in our subsequent analysis. Applying a standard torus 
model to any model of the luminosity function reproduces, 
by design, the mean offset between the B-band and hard X- 
ray luminosity functions, as the parameters of this model 
are tuned to reproduce this offset. As many observations 



show, the fraction of broad-line quasars increases with 


lumi- 


nosity (Steffen et al. 2003; Ueda et al. 2003; Hasinaet 


2004; 


ISazonov & Revnivtsev 2004; Baraeretal. 


2005; Simnson 


12005.) ■ and so reproducing the relationship between B-band 



and hard X-ray luminosity functions requires adding param- 
eters to the standard torus model which allow luminosity- 
dependent scalings, i.e. the class of "receding torus" models. 
These, again by construction, reproduce the distinction be- 
tween hard X-ray and B-band quasar luminosity functions, 
including the dependence of this difference on luminosity. 
These are, however, phenomenological models designed to 
fit these observations. Our simulations, on the other hand, 
provide a self-consistent description of the column density, 
which predicts the differences between hard X-ray, soft X- 
ray, and optical luminosity functions without the addition of 
tunable parameters or model features designed to reproduce 
these observations. 

Our fits are accurate down to low luminosities, as is clear 
from our prediction for the X-ray luminosity function at bolo- 
metric luminosities lO^L©. Furthermore, we have calcu- 
lated the predicted z < 0.1 luminosity function in the B-band 
as well as in Ha emission, using the conversion between the 
two from Hao et al. (2005) and comparing directly to their 
luminosity functions for Seyfert galaxies and low-luminosity 
active galactic nuclei (AGN) (both type I and II), and find 
that our distribution n(Lpeak) and model for quasar lifetimes 
and obscuration reproduces the complete observed luminosity 
function down to a B-band luminosity Mb ^ -16. Although 
our prediction falls below the observed Seyfert luminosity 
function at fainter magnitudes, there is no reason to believe 
that mergers should be responsible for all nuclear activity at 
these luminosities (and indeed alternative fueling mechanisms 
for such faint objects likely exist) - it is surprising, in fact, 
that this picture reproduces the observed AGN activity to such 
faint luminosities. 

Using the bolometric corrections of 'Elvis et al.l m- 
stead of Marconi et al. ( 2004) results in a significantly steeper 
cutoff in the luminosity function at high bolometric luminosi- 
ties, as the bolometric luminosity inferred for the brightest ob- 
served X -ray quasars is alm ost an order of magnitude smaller 
using the Elvis et al. corrections. However, this is be- 

cause the Elvis et al. (1994) bolometric corrections do not 
account for any dependence on luminosity, and further the 
quasars in the sam ple of Elvis et al. (1994) are X-ray bright 
dKlvis et alJl2'?OTl . whereas it has been well-established that 
the ratio of bolometric luminosity to hard or soft X-ray lumi- 
nosity increases_W2thJn2 Wilkes et M 
11^94; Green etal. 1995; Vignali et al. 2003; Strateva et af 
|2£)05). Recent comparisons between large samples of quasars 
selected by both optical and X-ray surveys ( Risaliti & Elvis 
l2005ii) further suggests that this is an intrinsic correlation, 
not driven by e.g. the dependence of obscuration on lumi- 
nosity. For a direct comparison of the bolometric luminos- 
ity functions resulting from the two corrections, we refer to 
iHopkins et al.l(l2005dh . Our analysis uses the form for the UV 



to X-ray flux ratio, aox, from lVignali et alJ (l2003h . but our re- 
sults are relatively insensitive to the different values found in 
the literature. It is important to account for this dependence, 
as it creates a significant difference in the high-luminosity end 
of the bolometric quasar luminosity function and implies that 
a non-negligible fraction of the brightest quasars are not seen 
in optical surveys (s ee the discussion in iMarconi et al ] l200i 
iRichards et alJ2005h . 

Finally, our fitted form for the evolution of the break lu- 
minosity, with cx expikiT), cannot continue to arbitrar- 
ily high redshift. At redshifts z > 2-3, this asymptotes be- 
cause r ^ 1, whereas the observed quasar population declines 
above z ^2. This difference is not important for most of our 
calculated observables, as they are either independent of high- 
redshift evolution or evolve with cosmic time in some fashion 
as cx Jn(Lpeak)df, with little time and thus negligible contri- 
butions to integrated totals at high redshifts. However, some 
quantities, in particular the high-mass end of the black hole 
mass function (see §|6j, which is dominated by the small num- 
ber of the brightest quasars at high redshifts, can receive large 
relative contributions from these terms. Therefore, it is impor- 
tant in estimating these quantities to be aware of the turnover 
in the quasar density at high redshifts. 

We quantify this in Figuref^ where we show the predicted 
broad-line luminosity function (where the broad-line phase 
is determined below in § |6} in six luminosity intervals from 
z ^ 1.2-4.8. The intervals are those o f the COMBO- 
17 luminosity function from IWol f et all ( l2003t) . but we 
further comp are to the observed luminosity functions of 
'Warren et al. (1994), 'Schmidt. Schneider. & Gunn (1995^ 
Kennefick, Diorgovski, & De Carvalho ( 1995), Fan et al] 
(2001), and Richards et al. (2005) at the appropriate (labeled) 
redshifts. At each redshift z > 2, we take the fitted n(Lpeak) 
distribution above ('Equations l2 1 1 122> and rescale it according 
to an exponential cutoff: either pure density evolution (PDF), 
«(Lpeak) n(ipeak) X 10""™*""^', or purc peak luminosity 
evolution (PPLE), U^Ux lO-^p^LEfc-a) pij^jng j^e data 
gives apDE ~ 0.65 and appLE ~ 0.55, {y^ /v « 1.3 for both) 
in reasonable agreement with the density evolution of e.g. 
Fan et al. (2001). We note that this evolution, extrapolated as 
far as z ^ 6, is consistent also with the constraints on z ~ 6 
quasars from Fa n et alJ (12003). especially in the PPLE case. 

In each panel, we plot the resulting broad-line luminos- 
ity function (see § |3, for both the minimum and maximum 
redshift of the redshift bin, and both the PPLE (solid lines) 
and PDF (dashed lines) cases. The degeneracy between these 
possibilities is well-known, as current observations do not 
resolve the break in the luminosity function. Furthermore, 
the predicted luminosity function should be considered un- 
certain especially at low luminosities, as the quasar lifetime 
at these luminosities and redshifts can become comparable to 
the age of the Universe, at which point our formalism for the 
luminosity function as a function of «(Lpeak) becomes inac- 
curate. However, we are able to make testable predictions, 
based on differences between the two models in integrated 
galaxy properties (for example, color-magnitude diagrams of 
red sequence galaxies at low masses or the fraction of recently 
formed spheroids as a function of mass and redshift), which 
distinguish the PPLE and PDF models for the evolution of th e 
quasar luminosity function at z > 2-3 ("Hopk ins et al.l2005el) . 
Owing to these degeneracies and the poor constraints on the 
observed high-redshift luminosity functions, we have not con- 
sidered them (those at z > 3) in our fits to /i(Lpeak), but use 
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Fig. 1 1. — Running our predicted broad-line luminosity function (determined in to high redshifts, with either total density (dashed lines) or break 

luminosity (L, ; solid lines) decreasing exponentially with redshift above z = 2. In each panel, our prediction is shown for the minimum and maximum 
redshift of the corresponding interval from the COMBO-17 luminosity function of Wolf et al. ( 2003) (W03; black squares). Other references for the obser- 
vations shown are: R05 - Richards et al. (2005), WHO - Warren etal. |1994), FOl - Fan et al. i 2001 ), SSG - Schmidt. Schn eider, & Gunn 1 1995). KDC - 
IKennefick. Diorgovski. & De Carvalhoitl995i) . 



them here to roughly constrain the turnover in the quasar den- 
sity above z ~ 2 (i.e. fitting to opde and opple)- Which form of 
the turnover we use makes little difference in our subsequent 
analysis, but, as discussed above, including some turnover is 
important in calculating select quantities such as the extreme 
high-mass end of the black hole mass function. 

3.3. The Observed Distribution 

Given the column density distributions and quasar life- 
times calculated from our simulations in §|2] and the quantity 
n(i'peak) determined above ('§ I3.2> . we can predict the distri- 
bution of column densities observed in a given sample. This 
will depend not only on the range of observed luminosities 
and the redshift of the sample, but also on the minimum ob- 
served magnitude and frequency (i.e. the selection function) 
of the sample. For a nearly complete sample or estimate 
of the luminosity function, for example the hard X-ray lu- 
minosity function, at least to A^h ~ 10^^ cm~^, we can inte- 
grate the A^H(i',i'peak) distribution over the n'(i'peak) distribu- 
tion (weighted by the lifetime at L). 

Figure [21 plots the resulting distribution of column densi- 
ties for this analysis. The left pan el reproduces an d expands 
upon a portion of Fig. 3 of Hopk ins et alJ (I2005bl) . showing 
the distribution of column densities (scaled linearly) expected 
from the characteristic quasars Lpeak ~ i-* of the luminosity 
function observed in optical samples, based on the simulated 
column density distributions as a function of luminosity and 
peak luminosity (solid black line). Specifically, we plot the 



distribution of neutral A^hi values requiring that the observed 
B-band luminosity be above some reference value Lb. mm- The 
smooth curve shown is the best-fit to the Eg-v distribution 
of bright SDSS quasars with z < 2.2, from Hopkins et all 
(2004). The curve has been rescaled in terms of the col- 
umn density (inverting our gas-to-dust prescription) and plot- 
ted about a peak (mode) A^hi (undetermined in Hopkins et 
al. 2004) of A^Hi ~ 0.5 x 10^' cm'l The observationally im- 
plied Eb-v distribution is determined from fitting to the dis- 
tribution of photometric reddening in all SDSS bands (i.e. 
using the five-band photometry as a proxy for spectral fit- 
ting) in Sloan quasars, relative to the modal quasar colors at 
each redshift, for quasars with an absolute magnitude limit 
Mi < -22. The /-band absolute magnitude limit imposed in 
the observed sample, M, < -22, corresponds approximately to 
our plotted B-band limit ob,s > 10" . This estimate does 
not account for bright but strongly reddened quasars having 
their colors altered to the point where color selection crite- 
ria of quasar surveys will not include them. However, this 
effect would only serve to bring our distribution into better 
agreement with observations, as it would slightly lower the 
high-A^Hi tail. We also consider the predictions of a standard 
torus model and receding (luminosity-dependent) torus model 
in the figure (dashed and dotted lines, respectively). These 
should not be taken literally in this case - they reflect that 
these phenomenological models do not predict the distribu- 
tion of low/moderate column densities, but rather assume that 
all lines of sight not intersecting the torus are "unobscured," 
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Fig. 12. — Left panel: Distribution of column densities expected from the characteristic quasars Lp^ak ~ i» of the luminosity function observed in optical 
samples, for a standard torus model of quasar obscuration (dashed), a receding torus model (dotted), and the distributions of column densities as a function 
of instantaneous and peak luminosity in our simulations (solid). The distribution of neutral A'hi values is obtained requiring an observed B-band luminosity 
> 10" L0. The smooth red curve is the best-fit to the Eb-v distribution of bright SDSS quasars with z < 2.2, from Hopkins et alj 120041) . rescaled to column 
densities and plotted about a peak (mode) Nwi (undetermined in Hopkins et al. 2004) of Nwi fsj 0.5 X 10^' cm"-. The r'-band absolute magnitude limit imposed in 
the observed sample, M, < -22, corresponds approximately to our plotted B-band limit Lg obs > 10" Lq . Reproduced from Hopkins et al. 1 2005b). Right panel: 
Integrated distribution of total (neutral and ionized) column densities expected for a complete hard X-ray sample, from the column densities of our simulations 
and the «(Lpeak) distribution. The distr ibution below 10^* c m~^ is shown (dot- dashed hne) and re-p lotted as a single bin at Nw = 10'" cm"- for our modeled 
columns. Data shown are the results of lTreister et alj i2004l) (blue squares) and lMainieri et alj 120051) (red circles), with assumed Poisson errors. Solid squares 
assume an intrinsic photon index F = 1.9, for the soft X-ray quasar spectrum, open squares F = 1.7. 



and encounter some constant, small column density (usually 
chosen to be A^h ^ 10^" cm"^). 

The right panel of shows the integrated distribution (in 
logA^n) for a complete hard X-ray sample, both as predicted 
from our simulations based on the joint distribution of col- 
umn density, luminosity, and peak luminosity (solid), and for 
both the standard torus mod el (d ashed) and receding torus 
model (dot ted) described i n § 12.51 The data shown are the re- 
sults of Tre ister et al. !(^2004) (blue squares) and Mainieri et ajj 
( 12005) (red circles), with assumed Poisson errors, from multi- 
band Chandra and HST observations of GOODS fields. The 
solid squares are obtained by assuming an intrinsic photon in- 
dex for the soft X-ray quasar spectrum of F = 1 .9, the open 
squares assuming F = 1.7. For the sake of direct comparison 
with observed distributions, objects with A^h < 10^' cm~^, for 
which only an upper limit to the column density would be 
determined in X-ray observations, are grouped together and 
plotted as a single bin at A^h = 10^° cm"^. The actual distri- 
bution below 10^' cm"^ is shown as a dot-dashed line. We 
note that our model of the quasar spectrum assumes a photon 
index F = 1 .9 in the soft X-ray, but this has no effect on the 
column densities calculated from the surrounding gas in our 
simulations. 

The agreement between the observed column density dis- 
tribution and the result of our simulations once the same se- 
lection effect is applied supports our model for quasar evo- 
lution, and the good agreement extends to both optical and 
X-ray samples. Probing to fainter luminosities or frequen- 
cies less affected by attenuation broadens the column den- 
sity distribution, as is seen from the inferred column den- 
sity distributions in the X-ray. This broadening occurs be- 
cause, at lower luminosities, observers will see both intrin- 



sically bright periods extinguished by larger column densi- 
ties (broadening the distribution to larger A^e values) and in- 
trinsically faint periods with small column densities (broad- 
ening the distribution to smaller A^h values). The distribu- 
tion as a function of reference luminosity is a natural conse- 
quence of the dynamics of the quasar activity. Throughout 
much of the duration of bright quasar activity, column den- 
sities rise to high levels as a result of the same process that 
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bright quasars strongly reddened by large A^hi- Furthermore, 
a significant number of quasars are extinguished from optical 
samples or attenuated to lower luminosities, giving rise to the 
distinction between luminosity functions in the hard X-ray, 

soft X-ray, and optical. 

The standard torus model described in § 12.51 although un- 
able to predict the distribution of column densities seen in op- 
tically, relatively unobscured quasars, does a fair job of repro- 
ducing the observed distribution of X-ray column densities. 
The parameters of the model are, of course, chosen to repro- 
duce the data sh own (the model parameters are taken from 
iTreister et a l.'2004). Nevertheless, our prediction is still abet- 
ter fit to the observed distribution, with /v k,2 as opposed 
to /v ~ 7 (although the absolute values depend on the es- 
timated systematic errors in the column density estimations). 
The receding torus model fares even more poorly in reproduc- 
ing the observed column density distributions, and is ruled out 
at high significance {y^ jv w 10), although this can be allevi- 
ated if the observed samples are assumed to be incomplete 
above A^H ^ lO^-'cm"^. This disagreement results because, in 
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order to match the observed scaling of broad-Une fraction with 
luminosity (see §|4]below), this model assumes a larger cover- 
ing fraction for the torus at lower luminosities, normalized to a 
similar obscured fraction as the standard torus model near the 
break in the observed quasar luminosity function. However, 
since quasars with luminosities below the break dominate the 
total number counts, this predicts that the cumulative column 
density distribution must be significantly more dominated by 
objects with large covering angles, giving a larger Compton- 
thick population, inconsistent with the actual observed col- 
umn density distribution. 

Although we do not see a significant fraction of extremely 
Compton-thick column densities A^h ^ 10^^ cm"^ in the dis- 
tributions from our simulations, our model does not rule out 
such values. It is possible that bright quasars in unusually 
massive galaxies or quasars in higher-redshift, compact galax- 
ies which we have not simulated may, during peak accre- 
tion periods, reach such values in their typical column den- 
sities. Moreover, as our model assumes ^ 90% of the mass of 
the densest gas is clumped into cold-phase molecular clouds, 
a small fraction of sightlines will pass through such clouds 
and measure column densities simila r to those shown for th e 
"cold phase gas" in, e.g. Figure 2 of iHopkins et'aP (HoOli), 
Nh > 102^-2''cm-2. 

Furthermore, we have not determined the "shape" at any 
instant of the obscuration (e.g. the dependence of obscura- 
tion on radial direction), as in practice, for most of the most 
strongly obscured phases in peak merger activity, the central 
regions of the merging galaxies are highly chaotic. Gener- 
ally, the scale of the obscuration in the peak merger phases 
is ~ lOOpc, quite different than that implied by most tra- 
ditional molecular torus models, but we note that our res- 
olution limits, ~ 20pc in the dense central regions of the 
merger, prevent our ruling out collapse of gas in the central 
regions into a smaller but more dense torus. However, sev- 
eral efforts to modjl_ traditional tori through radiative transfer 
simulations re.g.. iGranato & D anese 1994; Schar tmann et al.l 
I2OO5.) suggest significant column densities produced on scales 
of ~ 100-200pc, comparable to our predictions, and we note 
that only the solid angle covered by a torus, not the abso- 
lute torus scale, is constrained in the typical phenomenologi- 
cal torus model (e.g., Antonuccll993l) . 

Whether the obscuration of bright quasars originates on 
larger scales than is generally assumed is observationally 
testable, either through direct probes of polarized scat- 
tered hght tracing the obscuring/reflecting structure (e.g., 
IZakamska et 31.1 *2005). or through correlations between ob- 
scura tion and e.g. host g alaxy morphologies and inclinations 
(e.g., iDonlev et al.ll2005l) . These larger scales typical of the 
central regions of a galaxy are widely accepted as the scales of 
obscuration in starbursting systems (e.g. Soifer et al. 1984a,b; 
Sanders et al. 1986, 1988a,b; for a review, see e.g. Soifer et 
al. 1987), which in our modeling is associated with rapid ob- 
scured quasar growth and precedes the quasar phase. Thus, 
it is natural to associate obscuration with these large scales in 
any picture which associates starbursts and rapid black hole 
growth or quasar activity, as opposed to the smaller scales 
~ pc implied by torus models primarily developed to repro- 
duce observations of quiescent, low-luminosity Type II AGN, 
which are usually not directly associated with merger activity. 
These low-luminosity AGN are in a relaxed state, suggesting 
the possibility that the remaining cold gas in the central re- 
gions of our merger remnants will collapse once the violent 
effects of the merger and bright quasar phase have passed. 



producing a more traditional small torus in a quiescent nu- 
cleus. The central point is that regardless of the form of obscu- 
ration, the typical magnitude of the obscuration is a strongly 
evolving function of time, luminosity, and host system prop- 
erties, and the observed column density distributions reflect 
this evolution. 

4. BROAD-LINE QUASARS 

4.1. Determining the Broad-Line Phase 

Optical samples typically identify quasars through their col- 
ors, relying on the characteristic non-stellar power-law con- 
tinua of such objects. However, observations of X-ray se- 
lected AGN show a large population of so-called Type 2 
AGN, most of which have Seyfert-like luminosities and typ- 
ical spectra in X-ra ys and wavelengths longward of 1 /^m 
(e.g., Elvis et al. 1994), but are optically obscured to the point 
where no broad lines are visible. Their optical continua, 
in other words, resemble those of typical galaxies and thus 
they are not identified by conventional color selection tech- 
niq ues in optical qua sar surveys. Traditional unification mod- 
els jAntonuccfll993h have postulated a static torus as the ex- 
planation for the existence of the Type 2 population, with 
such objects viewed through the dusty torus and thus op- 
tically obscured. Moreover, both synthesis models of the 
X-ray backgrou nd (Setti & Woltier 1989; Madau et al. 1994 
"Comastriet al. 'l995; Gifli et al. 1999, 2001) and recent di- 
rect observations in large surveys (e.g., Zakamska et al. 2004 
2005 ) indicate the existence of a population of Type 2 quasars, 
with similar obscuration but intrinsic (unobscured) quasar- 
like luminosities. 

Observations of both radio-loud ( Hill, Goodrich, & DePov] 
1996; Simcson. Rawlings. & Lacy 1999; Willott et al. 200(1 
Simcson & Rawlings 2000; Grimes. Rawlings. & Willora 
2004.) and radio-quiet (Steffen et al. 2003; Ueda et al. 20031 
Hasinge [^004; S azonov & Revnivtsev 2004; Barge r et alJ 
2OO5F ISimpson. 2005) quasars, however, have shown that 
the broad-line fraction increases with luminosity, with 
broad-line objects representing a large fraction of all AGN 
at luminosities above the "break" in the luminosity function 
and rapidly falling off at luminosities below the break. 
Modifications to the standard torus unification model ex- 
plain^ this via a luminosity-dependent inner torus radius 
llLawre nce 1991), but this represents a tunable modification 
to a purely phenomenological model. Furthermore, as the 
observations have improved, it has become clear that even 
these luminosity-dependent torus models cannot produce 
acceptable fits to the broad line fraction as a function of 
luminosity (e.g.. lSimpsonll2005l) . However, we have shown 
above that the obscuring column, even at a given luminosity, 
is an evolutionary effect, dominated by different stages of 
gas inflow in different merging systems giving rise to varying 
typical column densities, rather than a single static structure. 
It is of interest, then, to calculate when quasars will be 
observed as broad-line objects, and to compare this with 
observations of broad line quasars and their population as a 
function of luminosity. 

Figure [O] shows the B-band luminosity as a function of 
time for both the quasars and host galaxies in three rep- 
resentative simulations^ the A2, A3, and A5 cases de- 
scribed in detail in § 12.11 These simulations each have 
/gas = 1-0, ^Eos = l-0,Zgai = 0, with virial velocities Vyir = 
113, 160 and 320 km s"', with resulting final black hole 
masses M^^ = 3 x 10^, 3 x 10^, and 2 x \QPMq, respec- 
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Fig. 13. — Intrinsic (right panels) and median attenuated (left panels) B-band luminosity of the quasar (thick line) and host galaxy (integrated ove r all stars, 
thin line; ignoring bulge stars, dotted line) as a function of time. Results are shown from three representative simulations: A2, A3, and A5 (see § 12.11 with 
lEOS = 1-0, Zgai = 0, and virial velocities Vyji = 113, 160, and 320 kms"' . Each quasar should be observable as a broad-line AGN when Lb.qso }t ^B.host- Colors 
show the stellar light curve with different gas fractions /gas = 1.0 (black), /gas = 0.4 (blue), and /gas = 0.2 (red); quasar light curves are similar for each gas 
fraction. 



tively. The thick line in each case shows the quasar B- 
band luminosity, and the thin line shows the integrated B- 
band luminosity of all stars in the galaxy. New stars are 
formed self-consistently in the simulations according to the 
ISM gas pro perties, equation of state and st ar formation model 
described in ISpringel & Hernauisll J200 3'). with the age and 
metallicity taken from the local star-forming ISM gas, which 
is enriched by supernova feedback from previous star forma- 
tion. We then use the stellar population synthesis model of 
Ifiruzual & Chariot (2003) to determine the B-band luminos- 
ity (the B-band mass-to-light ratio) of new stars based on the 
stellar age and metallicity. The dotted line shows the result 
neglecting bulge particles, which must be initialized at the 
beginning of the simulation with random or uniform ages and 
metallicities instead of those quantities being determined self- 
consistently from the simulation physics. The right panels 
plot the intrinsic values of these quantities, and the left panels 
plot the median observed values of these quantities, where we 
have used our met hod f or determining column densities and 
dust attenuation (§ 12. 2> to every star and bulge particle for 
each line of sight. 

Unfortunately, the host galaxy luminosity does not scale 
with instantaneous and peak quasar luminosity as do, for ex- 
ample, the quasar lifetime and obscuration. Rather, there are 
important systematic dependencies, the largest of which is the 
dependence on host galaxy gas fraction. If the host galaxies 
are more massive, more concentrated, or have a weaker ISM 
equation of state pressurization, then they will more effec- 



tively drive gas into the central regions and maintain high gas 
densities for longer periods of time, as the deeper potential 
well or lack of gas pressure requires more heat input from the 
quasar before the gas can be expelled. These conditions will 
generally produce a quasar with a larger peak luminosity (fi- 
nal black hole mass), but also form more new stars, meaning 
that the B-band relation between host and quasar luminosity 
is roughly preserved. 

However, the the black hole consumes only a small fraction 
of the available gas (comparison of e.g. the stellar mass and 
black hole mass suggests the black hole consumes ~ 0.1% of 
the gas mass), and so, at least above some threshold /gas ^0.1, 
the quasar peak luminosity does not significantly depend on 
the galaxy gas fraction (see, e.g. Figure 2 of Robertson et 
al. 2005b). But, the mass of new stars formed during the 
merger does strongly depend on the available gas. For exam- 
ple, simulations which are otherwise identical but have initial 
/gas = 0.2, 0.4, 0.8, 1.0 (i.e. an increasing fraction of the ini- 
tial disk mass in gas instead of stars) produce similar peak 
quasar luminosity and final total stellar mass (within ^ 30% 
of one another), reflecting the conversion of most gas into 
stars and the fact that the peak quasar luminosity is deter- 
mined more by the depth of the potential well than the total 
available gas supply. But, the mass of new stars formed in a 
merger scales roughly as M,j.new oc /gas (as it must if the initial 
gas fraction does not change the final total stellar mass), and 
since young stellar populations dominate the observed B-band 
luminosity (especially during the peak merger and starburst 
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Fig. 14. — Ratio of observed (attenuated) B-band quasar luminosity to host 
galaxy luminosity as a function of the ratio of instantaneous to peak quasar 
bolometric luminosity. Results are from simulations A2 (black diamo nds), 
A3 (blue circles), and A5 (red X 's) (the same simulations shown in Figure lTsl 
with qEOS = 10, Zsal = 0, and virial velocities Vyir = 113, 160, and 320 km s"' . 
Each panel shows the same simulations except for a different initial gas frac- 
tion /gas = 0.2, 0.4, 0.8, 1.0 as labeled. Solid lines are the predictions of 
Equation l26l 



phases associated with the bright quasar phase of interest), 
this implies roughly that Lb cx /gas- 

We demonstrate this explicitly in Figure[0] where we show 
in each panel the host galaxy and stellar B-band light curves 
for otherwise identical simulations with different gas frac- 
tions, /gas = 0.2 (red), 0.4 (blue), and 1.0 (black). In each 
of these cases, the quasar light curve is nearly identical (we 
show only the /gas =1.0 quasar lightcurve, for clarity, but the 
others are within ^ 30% of the curve shown at most times, 
with no systematic offset). 

In order for a quasar to be classified as a "broad-line" ob- 
ject, the optical spectrum must be visible and identified as 
such in the observed sample. This is clearly related to the ratio 
of quasar to host galaxy luminosity, but the threshold for clas- 
sification is not obvious. In an X-ray or IR-selected sample, 
optical follow-up should be able to disentangle host galaxy 
light and identify quasar broad-line spectra with fluxes a fac- 
tor of several fainter than the host. However, automated op- 
tical selection based on color or morphological criteria might 
well exclude objects unless the quasar luminosity is a factor 
of several greater than that of the host galaxy. Therefore, there 
is significant systematic uncertainty in the theoretical defini- 
tion of a broad-line quasar To first order, based on the above 
arguments, we can classify "broad-line quasars" as objects in 
which the quasar optical luminosity is larger than some mul- 
tiple /bl of the host galaxy optical luminosity. Because the 
relevant ratio is different depending on the survey and selec- 
tion techniques, we consider the range /bl = 0.3 - 3, with a 
rough median /bl = 1 . Furthermore, because our simulations 
do not allow us to model the broad-line regions of the quasar 
or spectral line structures as influenced by e.g. reddening and 
dust absorption, we adopt the B-band luminosity of the quasar 
and host galaxy as a proxy for optical luminosity and more 
complex (but often quite sample-specific) color and morpho- 
logical selection criteria. 

In Figure [O] the B-band host galaxy luminosity is quite 
flat as a function of time, relative to the quasar B-band lumi- 
nosity, and is roughly given by /Lq ^ M»,new/A^0, where 
A^*.new is the mass of new stars formed in the merger As 
noted above, this scales approximately Unearly with initial gas 



fraction at fixed final total stellar mass M», giving /Lq sa 
Cc,ai(M,/M0)/gas, where Cgai is a correction of order unity 
which we can fit from the simulations (essentially a mean 
mass-to-light ratio for the newly formed stars). The bolomet- 
ric correction of the quasar is usually defined by L^^" = cbL!^", 

and the quasar peak luminosity is Lpeak = cl LEdd(A^BH)' where 
again cl is a correction factor of order unity which we can cal- 
culate from our form for the quasar lifetime (see EauationllO> 
or measure in the simulations. 

If we require that the quasar B-band luminosity be larger 
than a factor /bl of the host galaxy B-band luminosity, we 
obtain 

LZ/Lq > /BLCBCgal(M,/MQ)/gas. (23) 

Dividing this through by Lpeak, we have 
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We can test this scaling relation against the results of our sim- 
ulations, and do so in Figure O Rearranging the equations 
above gives 
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which we can compare to our direct calculation of L'^° /L°g 
and Lboi /^peak for each simulation snapshot. 

Ultimately, we are not interested so much in the intrinsic 
B-band luminosity of the quasar and host galaxy, but rather 
the observed luminosities; i.e. we are interested in the ratio 
^TohJ^tohs = (ir/4'')(exp{-(Te-rc)}), where tq and 
Tc are "effective" optical depths which we use to denote the 
mean attenuation of quasar and host galaxy B-band luminosi- 
ties, respectively. We have considered the distribution of col- 
umn densities attenuating the quasar as a functi on o f instanta- 
neous and peak quasar luminosity in detail in § l2.3l above: the 
attenuation of the host galaxy as a function of luminosity, ob- 
served ban d, halo mass, and star formation rate are discussed 
in detail in lJonsson et alJ J2005h . Combining these fits gives, 

roughly, (exp{-(rQ-TG)}) ~ {M{j\(f MQf '^^ , but a better 
approximation can be determined directly from the simula- 
tions. 

This scaling can be understood roughly using toy models 
of uniformly mixed luminous sources within the galaxy de- 
scribed by J onsson et al. (.2005) , after accounting for the fact 
that the luminosity (star formation rate) dependent portion of 
the attenuation scales with luminosity in a similar manner 
to our quasar attenuation (compare our tq cx A^h oc Lq^o"'""'^'* 
to their (x Lg^^^j). The key consequence of this is that 
more massive systems (higher bulge and black hole masses) 
have their host galaxy light proportionally more attenuated 
in mergers, meaning that (as suggested by the comparison of 
light curves in Figure the quasar is more likely to be ob- 
served with an optical luminosity larger than that of its host. 

Figure O plots the ratio of the observed (attenuated) B- 
band quasar luminosity to the observed host galaxy B-band 
luminosity as a function of the ratio of instantaneous to peak 
quasar bolometric luminosity. We show the results for four 
different gas fractions /gas = 0.2, 0.4, 0.8, 1 .0 as labeled. For 
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each gas fraction, we consider our simulations A2 (black dia- 
monds), A3 (blue circles), and A5 (red x 's) (the same simula- 
tions shown in Figure[Oj with ^eos = 1-0, Zgai = 0, and virial 
velocities Vyir =113, 160 and 320kms~', using the labeled ini- 
tial gas fraction. The colored lines in each panel show the pre- 
dictions of combining the scalings expected for the intrinsic 
luminosities (Eauationl25> and attenuations as above, giving 
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where the colored lines each use the Mgjj and /gas of the simu- 
lation of the corresponding color and panel. This scaling pro- 
vides a good estimate of the observed optical quasar-to-galaxy 
luminosity ratio, including the complicated effects of atten- 
uation, evolving mass-to-light ratios, metallicities, and host 
galaxy properties, as a function of gas fraction, final black 
hole mass, and the ratio of the current to peak quasar lumi- 
nosity. Although, for clarity, we have not shown a range of 
simulations varying other parameters, we find that this scaling 
is robust to the large number of quantities we have considered 
in our simulations - there are systematic offsets in e.g. Lpeak 

andMgjj with changes such as e.g. different ISM equations of 

state, but the scaling in terms of Lpeak and Mg^ is unchanged. 

Because the ratio of observed quasar and host galaxy B- 
band luminosities in our simulations obeys the scaling of 
Equation|26l we can use it to predict the properties of "broad- 
line" quasars, defined by L^°^'^,^ > /blLI^'^^s- T° how- 
ever, we must assume a typical host galaxy gas fraction. Un- 
fortunately, because our empirical modeling in terms of the 
quasar lifetime as a function of L and Lpeak does not have a 
systematic dependence on host galaxy gas fraction (see § I2.4> . 
we have no constraint on this parameter. It is, however, con- 
venient for several reasons to consider /gas = 0.3 as a typical 
value for bright quasars. 

First, such a gas fraction is capable of yielding the brightest 
observed quasars; second, scaling a Milky-Way like disk with 
the observed z = gas fraction '--^ 0.1 to the re dshifts of peak 
quasar activity gives a similar gas fraction ('e.g.. lSpringel et aO 
12005 ah : third, gas fractions > 30% in major mergers are 
needed to explain the observed fundamental plane (Robert- 
son et al. 2005c, in preparation), kinematic properties (Cox 
et al. 2005c, in preparation), and central phase space den- 
sities (Hernquist, Spergel & Heyl 1993) of elliptical galax- 
ies; fourth, this choice implies that the brightest quasars with 
M^jj ^ 1O'°M0 attain observed B-band luminosities ^ 1000 
times that of their hosts at their peaks, as is observed (e.g., 
iMcLure & Dunlop 2004). Finally, and most important, the 
assumed /gas and /bl are degenerate in our predictions for the 
broad-line population, as they both enter linearly in the ra- 
tio of host galaxy to quasar B-band luminosity. Therefore, 
the range of /bl = 0.3-3 which we consider (for a fixed me- 
dian /gas = 0.3) can be equivalently considered, for a fixed 
median /bl = 1, to represent a theoretical uncertainty in the 
host galaxy gas fraction, /gas = 0.1-0.9; i.e. spanning the 
range from present, relatively gas-poor Milky-Way like disks 
to almost completely gaseous disks. This, then, gives for our 
"broad-line" criterion. 
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The "broad-line" phase is thus, as is clear from Figure [O] 
and implicit in our definition of the broad-Une phase, closely 



associated with the final "blowout" stages of quasar evolu- 
tion, when the mass of the quasar reaches that correspond- 
ing to its location on the Mbh - cr relation and gas is expelled 
from the central regio ns of the galaxy, shutting down accretion 
(Di Matteo et al. 2005). We note that combining the equation 
above with our fitted quasar lifetimes gives an integrated time 
when the quasar would be observable as a broad line object 
of fBL ~ 10-20Myr, in good agreement with the optically 
observable bright quasar lifetimes we calculate directly from 
our quasar light curves, including the effects of attenuation, 
and with empirical estimates of the quasar lifetime which are 
based directly on optically-selected, broad-line quasar sam- 
ples. 

The (Mgj]/1O^M0)"'^ term in the above equation reflects 
the fact that, below a certain peak luminosity, quasars are less 
likely to reach luminosities above that of the host galaxy, as 
can be seen in the uppermost panels of Figure for a final 
black hole mass of = 3x10^- i.e. the smallest AGN ai-e 
proportionally less optically luminous than their hosts. This 
does not imply that such systems are not inherently broad- 
line objects, but only that the host galaxy light will increas- 
ingly dominate at lower luminosities. We also caution against 
extrapolating this to large or small Mgj^, as the attenuation 
becomes more difficult to predict at these peak luminosities, 
and the linear formula above is not always accurate (see Fig- 
ureinj. 

We can use this estimate of the broad-line phase and our 
model of the quasar lifetime to calculate the total energy 
radiated in this bright, optically observable stage following 
the calculation of § 12.41 but with a minimum luminosity de- 
termined by Equation ^\ This gives an integrated fraction 
- 0.3-0.4 (- exp{-0.2/BL(/gas/0.3)/aL}) of the total radi- 
ant energy emitted during the broad-line phase. Thus, despite 
the short duration of this optical quasar stage, a large fraction 
of the total radiated energy is emitted (as it represents the fi- 
nal e-folding in the growth of the black hole) when most of 
the final black hole mass (§ 12. 4^ is accumulated. Account- 
ing for the luminosity dependence of our bolometric correc- 
tions (with the optical fraction of the quasar energy increas- 
ing with bolometric luminosity) as well as the small fraction 
of objects observable at lower luminosities (with larger typ- 
ical obscuring column densities) increases this fraction to as 
much as ^ 0.6-0.7 for bright quasars. Therefore, despite 
the fact that the duration of the optically observable broad- 
line quasar phase may be 1/10 that of the obscured quasar 
growth phase, the changing quasar luminosity over this period 
and non-trivial quasar lifetime as a function of luminosity im- 
plies only small corrections to counting arguments such as 
that of Soltan (1982), which rely on the total observed optical 
quasar flux density to estimate the relic supermassive black 
hole density. 

4.2. The Broad-Line Fraction as a Function of Luminosity 

By estimating the time that a quasar with some Lpeak will 
be observable as a broad-line quasar at a given luminosity, we 
can then calculate the broad-line quasar luminosity function 
in the sa me fa shion as the complete quasar luminosity func- 
tion in § 13.21 Instead of the full quasar lifetime df/dlogL, 
we consider only the time during which broad-lines would be 
observed (i.e. that the quasar spectrum would be recognized 
as opposed to the host galaxy spectrum), as identified in our 
simulations (§ 14. 1> . 

For a sample selected in hard X-rays (i.e. the selection func- 
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tion only being relevant at column densities > 10^"^ cm~^), we 
show the resulting "broad-line" luminosity function in Fig- 
ure|9](cyan dot-dashed lines), and compare it to the broad-line 
quasar luminosit y function ide ntified in the hard X-ray lumi- 
nosity function of lBarger et alJ ('2005). The agreement is good 
at all luminosities, and our model explains both the fact that 
broad-line quasars dominate the luminosity function at lumi- 
nosities well above the "break" in the luminosity function, and 
the downturn in the broad-line quasar population at luminosi- 
ties below the peak. Essentially, the broad-line quasar popula- 
tion more closely traces the shape of the n(Lpeak) distribution, 
giving rise to the observed behavior as a dual consequence 
of luminosity-dependent quasar lifetimes and the evolution- 
ary nature of quasar obscuration in our simulations. 

Figure [21 compares our theoretical predictions to the 
2dF-SDSS (2SLA0) g-band luminosity function of broad- 
line quasars from Richards et al. (2005) (black squares), as 
well as the B-band luminosity function from ICroom et al.l 
(|2004) (green circles), at several redshifts from z ~ 0.3 - 
2, over which range the surveys are expected to be rela- 
tively complete (for broad-line quasars). The 2dF-SDSS 
result is the most recent determination of the broad- 
line luminosity function, but compares well with previous 
determ in ations by, e.g., Bovle et al. (1988), Koo & Kron 
(11988"), Marano, Zamorani, & Zitelli (1988), Bovle et al. 
0990), Bovle. Jones. & Shanks (1991), ZitelH et al. (1992), 
iBovle et al. (2000), and Croometal. (2004). Open squares 
correspond to bins in luminosity which have been corrected 
for incompleteness following Page & Carrera (2000), but this 
correction is uncertain as the bins are not uniformly sampled. 
We compare this at each redshift to the prediction of our deter- 
mination of the quasar "broad-line" phase, where we estimate 
that the quasar is observable as a broad line object when its 
observed B-band luminosity is greater than a factor /bl = 1 of 
that of the host galaxy. We calculate this for both the mini- 
mum and maximum observed redshift of each bin to show the 
range owing to evolution of the luminosity function over each 
interval in redshift. The systematic uncertainty in our predic- 
tion can be estimated from the dotted lines, which show the 
prediction (at the mean redshift of the bin) if we instead re- 
quire the observed quasar B-band luminosity to be above a 
factor of 0.3 (upper lines) or 3 (lower lines) of the obser ved 
host galaxy B-band luminosity, which as discussed in § 14.11 
can alternatively be considered an uncertainty in host galaxy 
gas fraction, with /gas = 0.1 and /gas = 0.9, respectively. 

The agreement at all luminosities and redshifts shown is 
encouraging, given the simplicity of our determination of the 
broad-line phase from the simulations, but the systematic un- 
certainties are large, emphasizing the importance of calculat- 
ing detailed selection effects in contrasting e.g. "broad-line" 
samples from optical and X-ray surveys, as opposed to assum- 
ing a constant obscured fraction at a given luminosity based 
on the ratio of luminosity functions as has been adopted in 
previous phenomenological models. The difference between 
different choices of /gas is suppressed at the high luminos- 
ity (and correspondingly high redshift) end of the luminosity 
function, because the quasar-to-galaxy B-band luminosity ra- 
tio scales as oc (Mgjj)"'^; i.e. regardless of the choice of /bl, 
quasars increasingly overwhelm their host galaxy in large sys- 
tems near their peak luminosity. However, at low luminosity, 
the predictions rapidly diverge, implying that a measurement 
of the faint end of the broad-line quasar luminosity function, 
with a reliable calibration of /bl, can constrain the typical gas 



fractions of quasar host galaxies and the evolution of these gas 
fractions with redshift. 

By dividing out the predicted luminosity function 4)hx, we 
can estimate the fraction of "broad line" objects observed in 
reasonably complete X-ray samples as a function of luminos- 
ity. This is shown in Figure^] where for ease of comparison 
we have shown the broad-line fraction as a function of hard 
X-ray (2-10 keV) luminosity. Our prediction, based on de- 
termining the time a quasar with a given luminosity L and 
peak luminosity Lpeak in our simulations will be observable 
with a B-band luminosity greater than a fraction /bl = 1.0 
of the host galaxy observed B-band luminosity, is shown as 
th e thick black line. This is compared to the observations 
of lUedaetal l (l2003h (squar es), Hasinger (2004) (circles) 
'Grimes, Rawlines, & Willott (2004) (triangles), and Simpson] 
(2005) (diamonds). The data from Hasinger (2004) has been 
scaled from soft X-ray (0.5-2 keV) using our bolometric 
correc tions, and the data from iGrimes. Rawlines. & WillottI 
i2004l) andlSimpson (2005) have been converted from [O III] 
luminosity as in Simnson (2005) using the mean correction 
for Seyfert galaxies (Mulchaey et al. 1994), L[oiii] = 0.015 x 

i'2-lOkeV- 

We also plot as upper and lower dashed lines the results 
of changing /bl, the fraction of the host galaxy B-band lu- 
minosity above which the quasar B-band luminosity must be 
observed for identification as a "broad-line" object, consider- 
ing /bl = 0.3, and 3, respectively. We determine this for the 
low-redshift z < 0.3 quasar distribution, from which most of 
the data are drawn. The red dot-dashed line shows the differ- 
ence at high redshift, if just z > 1 quasars are considered (for 
/bl = !)■ The broad-line fraction is systematically lower, pri- 
marily because the break luminosity in the luminosity func- 
tion moves to higher luminosity with redshift, meaning that at 
a fixed luminosity below the break, a smaller fraction of ob- 
served objects are at L ^ ipeak in the "blowout" phase of peak 
optical quasar luminosity. Finally, the dotted line shows the 
results assuming a "light bulb" model for the broad-line phase 
(but still using our n(Lpeak) distribution, otherwise this trans- 
lates to a constant obscured fraction with luminosity) life- 
times, with a fixed broad-line lifetime of fg = 20Myr. 

The prediction of the most basic torus model, with con- 
stant broad-line fraction ~ 0.36, is ruled out to high sig- 
nificance ix^ /v = 18.5, 17.2 if we consider all data points, 
or if we consider only the most well-constrained data, from 
Simpson [2005], respectively). Furthermore, the solid cyan 
line shows the best-fit luminosity-dependent torus model, in 
which the broad line f raction is given by (e.g.. ISimpsonl 1 99*^ 
IGrimes. Rawlines. & WillotL2004) 

/=l-l/x/l+3L/Lo, (28) 

where Lq is the luminosity where the number of broad line ob- 
jects is equal to the number of non-broad line objects. This fit 
is at best marginally acceptable over a narrow range in lumi- 
nosities (y^ /i'= 14.0, 7.3). Modified luminosity-dependent, 
receding torus models have been proposed which give a better 
fit to the data by, for e xample, allow ing the torus height to vary 
with luminosity (e.g., Simpson 2005), but there is no physical 
motivation for these changes, and they introduce such varia- 
tion through additional free parameters that allow a curve of 
essentially arbitrary slope to be fitted to the data. 

However, the prediction of our model agrees reasonably 
well ix^ /v = 4.0, 1.2) with the observations over the entire 
range covered, a span of six orders of magnitude in luminos- 
ity. We emphasize that our prediction, which matches the 
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Fig. 15. — Broad-line quasar luminosity function oPRichards et al.'f2005^ from the 2dF-SDSS (2SLAQ) survey (black squares) and that o flCroom et al J 120041) 
(green circles) from the 2dF survey, compared to our predicted "broad-line" luminosity function from our determination of the relative quasar and host galaxy 
luminosities in our simulations (solid line), where we estimate that quasars are observable as "broad-line" objects when their observed B-band luminosity is 
greater than a factor /bl of that of the host galaxy. Solid lines are shown for the minimum and maximum observed redshift in each bin (as labeled), assuming 
/bl = 1 ■ Dotted lines show the result for /bl = 0.3 and /bl = 3, at the mean redshift of the bin, i.e. con'esponding to "broad-line" luminosity functions in surveys 
which are complete to quasars with observed optical luminosity ~ 0.3 and 3 times that of the host galaxy, respectively, or alternatively reflecting nearly complete 
theoretical uncertaint y regarding merging galaxy gas fractions (/gas = 0.1 —0.9). Open squares are observations with uncertain incompleteness corrections in 
IRichards etallOSna) . 



data better than standard torus models that are actually fit- 
ted to the data, is not a fit to the observations. Instead, it is 
derived from the physics of our simulations, including black 
hole accretion and feedback which are critical in driving the 
"blowout" phase which constitutes most of the time a quasar 
is visible as a "broad-line" object by our estimation, and from 
the n(Lpeak) distribution implied by our model of quasar life- 
times and the bolometric quasar luminosity function. The 
agreement suggests that our choice of the parameter combi- 
nation /fiL/gas = 0.3 is a good approximation. As noted above, 
this implies that calibrating /bl for an observed sample, com- 
bined with the mean broad-line fraction and our modeling, 
can provide a constraint (albeit model-dependent) on the host 
galaxy gas fraction of quasars at a given redshift, which can- 
not necessarily be directly measured even with difficult, de- 
tailed host galaxy probes, as gas is rapidly converted into stars 
throughout the merger The uncertainty plotted, while large, 
actually represents a larger theoretical uncertainty - as dis- 
cussed above, if an observational sample were well-defined 
such that it were complete to broad-line objects with observed 
optical luminosity above a fraction /bl of the host galaxy lu- 
minosity, the range we consider would correspond to a range 
/gas = 0.1-0.9 in the quasar host galaxy gas fraction, which 
the observations could then constrain. 

In our modeling, the broad line fraction as a function of 
luminosity does not depend sensitively on the observed lu- 



minosity function, as evidenced by the relatively similar pre- 
diction at high redshift. The evolution we do predict with 
redshift, in fact, agrees well with that found bv lBarger et alJ 
over the redshift range z = 0. 1 - 1 .2 (see also La Franca 
et al. 2005), an aspect of the observations which is not re- 
produced in any static or luminosity-dependent torus model 
but follows from the evolution of the quasar luminosity func- 
tion in our picture for quasar growth. However, we do cau- 
tion that gas fractions may systematically evolve with red- 
shift, and as discussed above, a higher gas fraction will give 
generally shorter "broad-line" lifetimes using our criteria of 
quasar optical luminosity being higher than some fraction of 
the host galaxy luminosity, which will also contribute to the 
evolution in the mean "broad-line" fraction with redshift. Fi- 
nally, neglecting the role of luminosity-dependent quasar life- 
times gives unacceptable fits to the data ix^ /v = 66.0, 77.5), 
as the broad-line fraction as a function of luminosity is a con- 
sequence of both the evolution of obscuration and the depen- 
dence of lifetime on luminosity. 

Our model for quasar evolution provides a direct physical 
motivation for the change in broad line fraction with luminos- 
ity and suggests that it is not a complicated selection effect. 
As an observational sample considers higher luminosities (i.e. 
approaches and passes the "break" in the observed luminos- 
ity function), a comparison of the luminosity function and the 
underlying n(Lpeak) shows that it is increasingly dominated by 
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Fig. 16. — Predicted "broad-line" fraction of a complete X-ray sample at 
low I < 0.3 redshift, from our simulations (where the object is observable 
as a broad-line quasar when it has an observed B-band luminosity greater 
than a factor /bl = 1.0 of that of its host galaxy), is shown (thick black 
line). The results, changing our /bl to 0.3 and 3.0 are shown, or equiv- 
alently of assuming a host galaxy gas fraction /ga.s = 0.1 or 0.9 instead of 
~ 0.3 (dashed), as are the results assuming a "hght bulb" model where 
quasars spend a fixed time tg = 20 Myr as broad line objects with a lumi- 
nosity of Lpeak (dotted). For comparison, the (scaled to 2-10 keV luminos- 
ity) observations of Ueda et al. |2003) (squares), Hasineei i2004) (circles), 
IGrim es. Rawlings. & WillotI i2004) (triangles), and Simpson i2005) (dia- 
monds) are shown. The predicted result at higher redshift (z > 1) is shown 
(red dot-dashed line), offset owing to the shift in break luminosity of the lumi- 
nosity function with redshift. The best-fit luminosity-dependent torus model, 
fitted to the data, is shown as the solid cyan line. The best-fit static torus 
model is a constant broad-line fraction ~ 0.3 (not shown for clarity). 




Fig. 17. — Predicted "obscured" fraction (solid line) in an X-ray sample 
with identical redshift and luminosity range to that of Ueda et al. 1 2003), as a 
function of hard X-ray (2-10 keV) luminosity. Vertical error bars show Pois- 
son errors estimated from the total time at a given luminosity across all our 
simulations (absolute values of the error bars should not be taken literally). 
The "obscured" fraction is defined as the fraction of quasars with X-ray col- 
umn densities A^ h > 10 ^^cm~~ in bins of AlogL2-iokcV- The observations 
from lUeda etal] l200l are shown as black squares. 



sources near their peak luminosity in the final stages of Ed- 
dington limited growth. The final stages of this growth ex- 
pel the large gas densities obscuring the quasar, rendering it a 
bright, optically observable broad-line object for a short time. 
Therefore, we expect that the fraction of broad-line objects 
should increase with luminosity in quasar samples, as indi- 
cated by the observations. 

Many observational measures do not consider a direct op- 
tical analysis of the quasar spectrum in estimating the frac- 
tion of broad-line objects as a function of luminosity. For 
example. Ueda et al. (2003) adopt a proxy, classifying as "ob- 
scured" any quasars with an X-ray identified column density 



A^H > 10^^ cm"^, and as "unobscured" quasars below this col- 
umn density. We can compare to their observations, using the 
column density distributions as a function of luminosity from 
our simulations, which cover the entire range in luminosity of 
the observed sample. Specifically, we use a Monte Carlo re- 
alization of these distributions, employing our fitted «(Lpeak) 
distribution at each redshift to produce a list of quasar peak 
luminosities and then generating all other properties based on 
the probability distribution of a given property in simulations 
with a similar peak luminosity. We describe this methodol- 
ogy in detail in § |8j and provide several such mock quasar 
distributions at different redshifts. 

In Figure we compare our estimated "obscured" and 
"unobscured" fractions as a function of hard X-ray luminos- 
ity, using the same definitions as well as redshift and lumi- 
nosity limits as the observed sample. The solid line shows 
our prediction, with vertical error bars representing Poisson 
errors, where the number of "counts" is proportional to the 
total time spent by simulations at the plotted luminosity (the 
absolute value of these errors should not be taken seriously). 
The "obscured" fraction is determined in bins of luminosity 
AlogL2-iokeV = 0.5. Despite our large number of simula- 
tions, there is still some artificial "noise" owing to incomplete 
coverage of the merger parameter space, namely the apparent 
oscillations in the obscured fraction. However, the mean trend 
agrees well with that observed (black squares), suggesting 
that the success of our modeling in reproducing the fraction of 
"broad line" objects as a function of luminosity is not a con- 
sequence of the definitions chosen above. We do not show the 
predictions of the standard and luminosity-dependent torus 
models, as (because essentially any line of sight through the 
torus encounters a column density A^h > 10^^ cm"^) the pre- 
dictions of these models are identical to those shown and com- 
pared to the same observations in FigurefTSl 

Our prediction that the fraction of broad-line objects should 
rise with increasing luminosity is counterintuitive, given our 
fitted column density distributions in which typical (median) 
column densities increase with increasing luminosity. This 
primarily owes to the simplicity of our A^h fits; we assume 
the distribution is lognormal at all times, but a detailed in- 
spection of the cumulative (time-integrated) column density 
distribution shows that at bright (near-peak) luminosities, the 
distribution is in fact bimodal (see e.g. Figure 3 of Hopkins 
et al. 2005b and Figure 2 of Hopkins et al. 2005d), represent- 
ing both the heavily obscured growth phase and the "blowout" 
phase we have identified here as the "broad line" phase. Over 
most of a simulation, we find the general trend shown in Fig- 
ure |3] and discussed above, namely that typical column den- 
sities increase with intrinsic (unobscured) luminosity. This 
is because the total time at moderate to large luminosities is 
dominated by black holes growing in the obscured/starburst 
stages; here, the same gas inflows fueling black hole growth 
also give rise to large column densities and starbursts which 
obscure the black hole activity. However, when the quasar 
nears its final, peak luminosity, there is a rapid "blowout" 
phase as feedback from the growing accretion heats the sur- 
rounding gas, driving a strong wind and eventually terminat- 
ing rapid accretion, leaving a remnant with a black hole sat- 
isfying the Mbh-ct relation. This can be identified with the 
traditional bright optical quasar phase, as the final stage of 
black hole growth with a rapidly declining density (allowing 
the quasar to be observed in optical samples), giving typi- 
cal luminosities, column densities, and lifetimes of optical 
quasars. In these stages, larger luminosities imply more vi- 
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olent "blowout" events, i.e. a brighter peak luminosity quasar 
more effectively expels the nearby gas and dust, rendering 
a dramatic decrease in column density at these bright stages 
(see Hopkins et al. 2005f). 

We are essentially modeling this bimodality in more detail 
by directly determining the "broad-line" phase from our sim- 
ulations. However, the broad line fraction-luminosity rela- 
tion we predict is also a consequence of the more compli- 
cated relationship between column density, peak luminosity, 
and bolometric and observed luminosity, as opposed to the 
predictions from a model with correlation between A^h and 
only observed luminosity. The key point is that we find, near 
the peak luminosity of the quasar, as feedback drives away 
gas and slows down accretion, the typical column densities 
fall rapidly with luminosity in a manner similar to that ob- 
served. In our model for the luminosity function, quasars be- 
low the observed "break" are either accreting efficiently in 
early stages of growth or are in sub-Eddington phases com- 
ing into or out of their peak quasar activity. Around and 
above the break, the luminosity function becomes dominated 
by sources at high Eddington ratio at or near their peak lumi- 
nosities. Based on the above calculation, we then expect what 
is observed, that in this range of luminosities, the fraction 
of objects observed with large column densities will rapidly 
decrease with luminosity as the observed sample is increas- 
ingly dominated by sources at their peak luminosities in this 
blowout phase. This also further emphasizes that the evolu- 
tion of quasars dominates over static geometrical effects in 
determining the observed column density distribution at any 
given luminosity. 

Finally, if host galaxy contamination were not a factor, we 
would expect from our column density model that, at low lu- 
minosities (L < IO'^Lq, well below the range of most obser- 
vations shown), the broad-line fraction would again increase 
(i.e. the obscured fraction would decrease), as the lack of gas 
to power significant accretion would also imply a lack of gas 
to produce obscuring columns. However, at these luminosi- 
ties, typical of faint Seyfert galaxies or LINERs, our model- 
ing becomes uncertain; it is quite possible, as discussed pre- 
viously, that cold gas remaining in relaxed systems could col- 
lapse to form a traditional dense molecular torus on scales 
~pc, well below our resolution limits. Furthermore, host 
galaxy light is likely to overwhelm any AGN broad-line con- 
tribution, and selection effects will also become significant at 
these luminosities. 

4.3. The Distribution of Active Broad-Line Quasar Masses 

Our determination of the "broad-line" or optical phase 
in quasar evolution allows us to make a further prediction, 
namely the mass distribution of currently active broad-line 
quasars. At some redshift, the total number density of ob- 
served, currently active broad-line quasars with a given Lpeak 
will be (in the absence of selection effects) 

nBL(ipeak) ~ "(i'peak) fBL(ipeak) , (29) 

where fBL(i'peak) is the total integrated time that a quasar with 
peak luminosity Lpeak spends as a "broad-line" object (us- 
ing our criterion for the ratio of the observed quasar B-band 
luminosity to that of th e host galaxy), given by integrating 
our formulae in § 14.11 or directly calculated from the sim- 
ulations. Since we have determined roughly that a quasar 
should be observable as a "broad-line" object at times with 
L > 0.2Lpeak primarily just after it reaches its peak luminos- 
ity, in the "blowout" phase of its evolution, we expect the in- 
stantaneous black hole mass at the time of observation as a 



broad-line quasar to be, on average, Mgjj w A^BH^^peak), where 

^BH ^ MEd d(Lpe ak) modulo the order unity corrections de- 
scribed in § 12.41 Using our fitted n(Lpeak) distribution from 
the luminosity function, extrapolated to low redshift (z ^ 0), 
and combining it with the integrated "broad-line" lifetimes 
from our simulations as above, we obtain the differential num- 
ber density of sources in a logarithmic interval in Lpeak- Fi- 
nally, we use our Equation[ro|forMgj^(Lpeak) determined from 
our fitted quasar lifetimes (demanding that Lrad = ^rMgj^c^) to 
convert this to a distribution in black hole mass. 

Our predicted n(MBH), i-e. the number of observed active 
quasars at low redshift in a logarithmic interval of black hole 
mass, is shown in Figure^] We consider the complete dis- 
tribution of active quasar masses, for both broad-line and non 
broad-line objects, in the left panel of the figure, and the dis- 
tribution of broad-line objects only, n{M^\(), in the right panel. 
On the left, we show the complete distribution which would 
be observed without any observational limits (dashed line). 
We calculate this from the distributions of Eddington ratios in 
our simulations, as a function of current and peak luminosity, 
and our fit to n(Lpeak) (as, e.g. for our Monte Carlo realiza- 
tions). We also consider the observed distribution if we apply 
the l uminosity limit for co mpleteness from the SDSS sam- 
ple of Heckman et alJJ200 4) (dotted), Lpiuj > 10'' L©, which 
using their bolometric corrections yields L > 3.5 x lO^L©, 
and then additionally applying their minimum velocity dis- 
persion cr > 70kms"' (dot-dashed). Finally, we can weight 
this distribution by luminosity (solid line) to compare directly 
to that determined in their Fig. 1. The red points are taken 
fr om the lu minosity-weighted black hole mass function of 
iHeckman et alJ |2004), which serves as a rough estimate of 
the active black hole mass distribution given their selection 
effects. Vertical error bars represen t the range in parame - 
terizations of the mass function from lHeckman et aiTj2004l) . 
including whether or not star formation is corrected for and 
limiting the sample to luminosities L > lO'^L© or Eddington 
ratios > 0.01. Horizontal errors represent an uncertainty of 
0.2 dex in the black hole mass estimation (representative of 
uncertainties in the Mbh - f relation used). The agreement is 
good, especially given the significant effects of the selection 
criteria and luminosity- weighting. 

We also consider the predictions of a "light-bulb" or "expo- 
nential / fixed Eddington ratio" model of the quasar lifetime 
for the active black hole mass distribution (red lines). For 
purposes of the active black hole mass function, the two pre- 
dictions are identical and independent of the assumed quasar 
lifetime (modulo the arbitrary normalization), as both assume 
that all observed quasars are accreting at a fixed Eddington 
ratio, giving the distribution of active black hole masses. The 
dashed line shows the prediction for the complete active black 
hole mass function, which rises sharply to lower luminosi- 
ties, as it must given a luminosity function which increases 
monotonically to lower luminosities. The solid line shows 
the prediction of such a model with the complete set of se- 
lection effects from Heckmanetal. (2004) described above 
applied, as with the solid black line showing the prediction 
of our modeling. Here, we chose the characteristic Eddington 
ratio « 1 .0 by fitting the predicted curve to the Heckman et afl 
(12004) observations. Note that both the characteristic Edding- 
ton ratio and lifetime (normalization) of the curve are fitted, 
so the relative normalization of this curve and our full model 
prediction are not the same; for example, the predicted to- 
tal absolute number of active Mbh > 10^ quasars is higher in 
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Fig. 18. — Predicted distiibution of currently active black hole masses, both considering all types (Type I & 11; left) and only those visible as broad-line 
quasars (Type I; right), at low z ^ 0.3 redshift, from our «(Lpeak) distribution and the estimation of the "broad-hne" phase directly from the simulations. In the 
left panel (all quasar types), we consider the result with arbitrarily faint luminosity Hmits (dashed line), and with the luminosity completeness limit (dotted) and 
both luminosity limit and velocity dispersion limit (dash-dot) of the SDSS sample of Heckman et al. ( 2004). We then consider the mass distribution with these 
limits, weighted by Olll luminosity, for direct comparison to the mass function of Heckman et al. 1 2004), shown as red circles (vertical errors represent the range 
in different parameterizations of the luminosity-weighted mass function from Heckman et al. ( 2004), their Fig. 1, horizontal errors a ~ 0.2 dex uncertainty in the 
black hole mass). Black lines show this for our full model, red lines show the full distribution (dashed) and distribution with the same weighting and selection 
effects as Heckman et al. ( 2004) (soHd) for a light-bulb or exponential light curve model of quasar evolution. At right, the distribution of active "broad-line" 
quasar masses (solid, where an object is a "broad-line" quasar if the observed quasar B-band luminosity is above a factor /bl = 1 of that of the host galaxy - 
dotted and dashed lines show the result if /bl = 0.3 or 3, respectively). Black lines show the prediction of the full model, red and blue lines the predictions of a 
light-bulb/exponential light curve model with a standard torus model (red) and receding torus model (blue) used to determine the broad-line fraction. 



the full model than in the light-bulb or exponential models. 
Still, it is clear that these models produce too broad a distri- 
bution of active black hole masses, in disagreement with the 
observations. We could, of course, obtain an arbitrarily close 
agreement with the observations if we fit to the distribution of 
accretion rates, but such a model would recover a quasar life- 
time and accretion rate distribution quite similar to ours, as 
is evident from the agreement between the predictions of our 
simulations and the observations. A purely empi rical model 
of this type is considered by e.g. ' Merlonil J2004ft . who finds 
that similar qualitative evolution in the quasar lifetime and 
anti-hierarchical black hole assembly to that predicted by our 
modeling is implied by the combination of quasar luminosity 
functions and the black hole mass function. 

On the right of the figure, we show our predicted mass dis- 
tribution for low-redshift, active "broad-line" quasars (solid 
black lines), where we estimate that an object is a "broad-line" 
quasar if the observed quasar B-band luminosity is above a 
factor /bl = 1 of that of the host galaxy - dotted and dashed 
lines show the result if /bl = 0.3 or 3, respectively, parameter- 
izing the range of different observed samples. As discussed 
above, the range of /bl shown can be, alternatively, thought 
of as a parameterization of uncertainty in the host galaxy gas 
fraction, if (in an observed sample), the sensitivity to see- 
ing quasar broad lines against host galaxy contamination is 
known. Therefore, the location of the peak in the active broad- 
line black hole mass function can be used, just as the mean 
broad line fraction vs. luminosity, as a test of the typical gas 
fractions of bright quasar host galaxies, and can constrain po- 
tential evolution in these gas fractions with redshift. 

The prediction shown is testable, but appears to be in 
good agreement with preliminary results for the distribution 



of active broa d-line black hole masses from the SDSS (e.g., 
iMcLure & D unlop 2004). The observations may show fewer 
low-mass black holes than we predict, but this is expected, as 
observed samples are likely incomplete at the low luminosi- 
ties of these objects (even at the Eddington limit, a lO^M© 
black hole has magnitude Mg ^ -16). If, in our model, we 
were to consider instead a standard torus scenario for the def- 
inition of the broad-line phase, we would predict the same 
curve as that shown in the left half of the figure (black dashed; 
our prediction for the cumulative active black hole mass func- 
tion). This is because the standard torus model predicts that 
a constant fraction of objects are broad-line quasars, regard- 
less of mass or luminosity, thus giving identical distribu- 
tions of Type I and Type II quasar masses. If we consider a 
luminosity-dependent or receding torus model, the prediction 
is nearly identical to the black line shown. This is because, as 
shown in Figure our prediction for the broad line fraction 
as a function of luminosity is similar to that of the receding 
torus model. The differences in the model predictions for the 
broad-line fraction as a function of luminosity do manifest in 
the prediction for the active broad-line black hole mass func- 
tion, but the difference in these models is smaller than the 
^ la range from different values of /bl shown. However, 
if we consider different models for the quasar light curve or 
lifetime, the predicted active broad-line mass function is quite 
different (as is the cumulative active black hole mass func- 
tion). 

We show the predictions of a light-bulb or exponential light 
curve model for quasar evolution in the figure, adopting either 
a standard torus model (red) or receding torus model (blue) to 
determine the broad-line fraction as a function of luminosity. 
For the standard torus model, this predicts that the broad line 
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mass function should trace the observed luminosity function, 
rising monotonically to lower black hole masses in power-law 
fashion (just as seen in the red dashed line in the left half of 
the figure for the cumulative black hole mass function). For 
the receding torus model, the active black hole mass func- 
tion shows a peak (because, at lower luminosities, there are 
more observed quasars, but a larger fraction of them are ob- 
scured). However, the location of this peak is at roughly an 
order of magnitude smaller black hole mass than for our pre- 
diction. This assumes a typical Eddington ratio ^ 1, which 
we have fitted to the cumulative black hole mass function - 
the peak in the broad-line active black hole mass function in 
these models could be shifted to larger black hole masses by 
assuming a smaller typical Eddington ratio, but this would 
only worsen the agreement with the cumulative black hole 
mass function of Heckman et al. (2004). Furthermore, a ro- 
bust difference between the models is that the light bulb or 
Eddington-limited/exponential models predict, for the stan- 
dard torus case, no turnover in the active broad-line black hole 
mass function, and for the receding torus case, a broader dis- 
tribution in active broad-line quasar black hole masses than 
is predicted in our modeling. Roughly, the lognormal width 
of this distribution in our model is ^ 0.6 dex, whereas the 
light-bulb or exponential light curve models have a distribu- 
tion with width ^1.0 dex. As noted above, we obtain a similar 
prediction if we adopt our full obscuration model instead of 
the receding torus model here. A determination of the range 
of active, broad-line quasar masses can, therefore, constrain 
quasar lifetimes and light curves. 

Our model makes an accurate prediction for the distribution 
of active black hole masses, even at z ^ where our extrapola- 
tion of the luminosity function is uncertain. It is important to 
distinguish this from the predicted relic black hole mass dis- 
tribution, derived in § |6] which must account for all quasars, 
i.e. n(Lpeak) integrated over redshift. We additionally find for 
broad-line quasars, as we expect from our prediction of the 
broad-line phase, that these objects are primarily radiating at 
large Eddington ratios, / ^ 0.2- 1, but we address this in more 
detail in §|5] The success of this prediction serves not only to 
support our model, but also implies that we can extrapolate to 
fairly low luminosities, even bright Seyfert systems at z ^ 0. 
This suggests that many of these systems, at least at the bright 
end, may be related to our assumed quasar evolution model, 
fueled by similar mechanisms and either exhibiting weak in- 
teractions among galaxies or relaxing from an earlier, brighter 
stage in their evolution. As we speculate in § |8j our descrip- 
tion of self-regulated black hole growth may also be relevant 
to fainter Seyferts, even those that reside in apparently undis- 
turbed galaxies. 

5. THE DISTRIBUTION OF EDDINGTON RATIOS 

In traditional models of quasar lifetimes and light curves, 
the Eddington ratio, I = L/L^m is generally assumed to be 
constant. Even complex models of the quasar population 
which allow for a wide range of Eddington ratios according to 
some probability distribution P(l) implicitly associate a fixed 
Eddington ratio with each individual quasar, and do not al- 
low for P{1) to depend on instantaneous luminosity or host 
system properties. However, this is a misleading assumption 
in the context of our model, as the Eddington ratio varies in 
a complicated manner over most of the quasar light curve. 
Furthermore, the integrated time at a given Eddington ratio 
is different in different systems, with more massive, higher 
peak luminosity systems spending more time at large (/ ~ 1) 
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Fig. 19. — Distribution of Eddington ratios (left panels) and instantaneous 
black hole mass (right panels) as a function of quasar bolometric luminosity 
for our fiducial Milky Way-like A3 simulation, with = 160km s"' and 
^pcak ~ 5 X IO'^Lq. The trend of an increasingly narrow Eddington ratio and 
mass distribution (concentrated at higher values) with increasing luminosity 
is clear. The result of applying an ADAF-type radiative efficiency coiTection 
at low accretion rates is shown (dashed) as well as the result of considering 
only times after the final merger, with Mbh ~ ^bh (dotted). 



Eddington ratios. 

The probability of being at a given Eddington ratio should 
properly be thought of as a conditional joint distribution 
P(l \ L. Lpeak) in both instantaneous and peak luminosity, just 
as the quasar "lifetime" is more properly a conditional dis- 
tribution f()(L|Lpeak)- Rather than adopting a uniform Ed- 
dington ratio or Eddington ratio distribution, empirical esti- 
mates must consider mor e detailed formulations such as the 
framework presented in Steed & Weinberg! J2003h . which al- 
lows for a conditional bivariate Eddington ratio distribution 
and can therefore incorporate these physically motivated de- 
pendencies and complications in de-convolving observations 
of the quasar luminosity function to determine e.g. Eddington 
ratio distributions, active black hole mass functions, and other 
physical quantities. 

Figure ^] shows the distribution of Eddington ratios as a 
function of luminosity for the fiducial. Milky Way-like A3 
simulation (Vvir = 160kms"'). Over the course of the simula- 
tion, the system spends a roughly comparable amount of time 
at a wide range of Eddington ratios from / ^ 0.001 - 1 . At high 
luminosities, L > lO'^L© for a system with Lpeak ~ lO'^'L©, 
the range of Eddington ratios, is concentrated at high values 
I ^ 0.5- 1 with some time spent at ratios as low as / ^0.1. 
Note, however, that the y-axis of the plot is scaled logarithmi- 
cally, so the time spent at Z ~ 0. 1 in this luminosity interval is a 
factor ^ 5 smaller than the time spent at / > 0.5. Considering 
lower luminosities 10" < L < lO'^L©, the distribution of 
Eddington ratios broadens down to I ^ 0.01. Going to lower 
luminosities still, L < 10" L©, the distribution broadens fur- 
ther, with comparable time spent at ratios as low as / ~ 0.001, 
and becomes somewhat bimodal. At large luminosities near 
Lpeak, the system is primarily in Eddington-limited or near- 
Eddington growth. However, as we consider lower luminosi- 
ties, we include both early times when the black hole is grow- 
ing efficiently (high I) and late or intermediate times when the 
black hole is more massive but the accretion rate falls (low 
I). As we go to lower luminosities, the total time spent in 
sub-Eddington states increasingly dominates the time spent at 
Z ^ 1 , although the time spent at any given value of I is fairly 
flat with log(/). 
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Roughly, at some luminosity L, there is a constant proba- 
bility of being in some logarithmic interval in I, 
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and P(/|L,Lpeak) = otherwise. This is especially clear if we 
compare the distribution of Eddington ratios in each luminos- 
ity range obtained if we consider only times after the final 
merger of the black holes (dotted histograms). At the high- 
est luminosities, the distribution is identical to that obtained 
previously, since all the time at these luminosities is during 
the final merger However, as we move to lower luminosi- 
ties, the characteristic I move systematically lower, as we are 
seeing only the relaxation after the final "blowout" near Lpeak, 
with characteristic Eddington ratio / = L/Lpeak at any given 
luminosity L. These trends are also clear if we consider the 
distribution of instantaneous black hole masses in each lumi- 
nosity interval shown in the figure, which is trivially related 
to the Eddington ratio distribution at a given luminosity L as 
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Of course, it is clear here that Mbh ~ = 3 x 10 if we 
consider only times after the final merger. 

It has also been argued from observations of stellar black 
hole binaries that a transition between accretion states oc- 
curs at a critical Eddington ratio m = M/MEdd, from ra- 
diativel y inefficient accretion flows at l ow accretion rates 
(e.g., Esin. McClintock. & Naravan 1997) to radiatively effi- 
cient accretion through a standard Shakura & Sunvaev ( 1973) 
disk. Although the critical Eddington ratio for supermassive 
black holes is unc ertain, observations of black hole binaries 
(|Maccaron j l20()l as well as theoretical extensions of ac- 
cretion models (e.g.. lMever. Liu. & Mever-H ofmeister''2000') 
suggest nhait ^ 0.01. We can examine whether this has a 
large impact on our predictions for the luminosity function 
and «(Lpeak) distribution, by determining whether the distri- 
bution of Eddington ratios is significantly changed by such a 
correction. Because we assume a constant radiative efficiency 
L = ErMc^ with e,- = 0.1, we account for this effect by multi- 
plying the simulation luminosity at all times by an additional 
"efficiency factor" /eff which depends on the Eddington ratio 
I = L/ LEdd, 

/ 1 if Z > 0.01 
J'^ft-]^ 100/ if / < 0.01. ' 

Thi s choice for the efficiency factor follows from ADAF mod- 
els jNaravan&Yil 11995) and ensures that the radiative effi- 
ciency is continuous at the critical Eddington ratio l„h = 0.01. 
Applying this correction and then examining the distribution 
of Eddington ratios as a function of luminosity (dashed his- 
tograms in Figure I19> . we see that the distribution of Ed- 
dington ratios is essentially identical, with only a slightly 
higher probability of observing extremely low Eddington ra- 
tios / < 0.001 . Of course, our modeling of accretion processes 
does not allow us to accurately describe ADAF-like accretion 
at these low Eddington ratios, but such low values are not rele- 
vant for the observed luminosity functions and quantities with 
which we make our comparisons. This implies that such a 
transition in the radiative efficiency with accretion rate should 
not alter our conclusions regarding the luminosity function 
and the n(Lpeak) distribution (essentially, the corrections are 
important only at luminosities well below those relevant in 
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Fig. 20. — Predicted distribution of Eddington ratios based on the lumi- 
nosity function and tlie quasar evolution in our simulations, in two redshift 
intervals z < 0.5 (upper left) and 1.5 < z < 3.5 (upper right). The observed 
distributions for radio lou d (black squares ) and radio quiet (green circles) 
quasars are shown from Vestergaard i'2()04') with Poisson errors. Thick black 
lines show the predicted distribution given the same minimum observed lumi- 
nosity as the observed sample. Thin red lines show the predicted distributions 
for a sample extending to arbitrarily faint luminosities, dotted lines show the 
same, with the ADAF correction of §|5]applied at low accretion rates. Blue 
dashed lines show the prediction for a fixed (luminosity-independent) Ed- 
dington ratio distribution in a light-bulb or exponential light curve model, 
fitted to the z < 0.5 data and used to predict the 1 .5 < z < 3.5 Eddington ratio 
distribution given the observational luminosity limit. Lower panels show the 
predicted distributions for z ^ 1 in two luminosity intervals, above and be- 
low the "break" luminosity in the observed luminosity function (red lines here 
correspond to an observed (attenuated) B-band luminosity Lg obs > 10" Z.q). 



constructing the observed luminosity functions; see also Hop- 
kins et al. 2005c for a calculation of the effects of such a cor- 
rection on the fitted quasar lifetime and n(Lpeak) distributions, 
which leads to the same conclusion). 

Despite the broad range of Eddington ratios in the simu- 
lations, this entire distribution is unlikely to be observable in 
many samples. The effect of this can be predicted based on the 
behavior seen in Figure^] For example, we consider the dis- 
tribution of Eddington ratios that would be observed if the B- 
band luminosity Lg.obs > 10" L©, comparable to the selection 
limits at high redshift of many optical quasar samples. As ex- 
pected from the change in / with luminosity, this restricts the 
observed range of Eddington ratios to large values / ^ 0. 1 - 1, 
in good agreement with the range of Eddington ratios actu- 
ally observed in such samples. Essentially, it has reduced the 
observed range to a bolometric luminosity L > lO'^L© in the 
case shown, giving a similar distribution to that seen in the 
lower panel of the figure. 

We compare our predicted distribution of Eddington ratios 
to observations in Figure|20| Using the distribution of peak lu- 
minosities n(Lpeak) determined from the luminosity function, 
we can integrate over all luminosities to infer the observed 
Eddington ratio distribution, 

P(/)cx J dlogL y^dlogLpeak 

XP(/|L,Lpeak)^^^^^«(ipeak)- (33) 

As our estimate of P(/|L,Lpeak) above is rough, we 
do this by binning in Lpeak and averaging the binned 
P(/|L,Lpeak)dr/dlogL for each simulation in the range of 
Lpsik, then weighting by n(Lpeak) and integrating. We con- 
sider both the entire distribution that would be observed in the 
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absence of selection effects (red histograms), and the distri- 
bution observed demanding a B-band luminosity above some 
reference value, L^ obs > L^in (black histograms). The results 
are shown for redshifts z < 0.5 and z = 1.5- 3.5 , along with 
the observed distribution from Vestergaard (2004), with as- 
sumed Poisson errors. The observations should be compared 
to the black histograms, which have luminosity thresholds 
L= 10"'Lq and 10" L© for z< 0.5 and z = 1.5-3.5, respec- 
tively, corresponding approximately to the minimum observ- 
able luminosities in the observed samples in each redshift in- 
terval. 

The agreement is good, given the observational uncertain- 
ties, and it suggests that the observed Eddington ratio distribu- 
tion can be related to the non-trivial nature of quasar lifetimes 
and light curves we model, rather than some arbitrary distribu- 
tion of fixed / across sources. However, the selection effects 
in the observed samples are quite significant - the complete 
distribution of Eddington ratios is similar in both samples, 
implying that the difference in the observed Eddington ratio 
distribution is primarily a consequence of the higher luminos- 
ity limit in the observed samples - and a more detailed test of 
this prediction requires fainter samples. 

Still, there is a systematic offset in the observed samples 
at z < 0.5 and z= 1.5-3.5 which does not owe to selection 
effects. At progressively lower redshifts, more quasars with 
luminosities further below the "break" in the luminosity func- 
tion are observed, and therefore the observed Eddington ratio 
is broadened to lower Eddington ratios I ^ 0.1, whereas at 
high redshift the distribution is more peaked at slightly higher 
Eddington ratios. This difference, although not dramatic, is a 
prediction of our model not captured in "light bulb" or "fixed 
Eddington ratio" models, even when allowing for a distribu- 
tion of Eddington ratios, if such a distribution is static. We 
demonstrate this by fitting the low-redshift Eddington ratio 
distribution to a Gaussian (blue dashed lines in upper left), 
and then assuming that this distribution of accretion rates is 
unchanged with redshift, giving (after applying the same se- 
lection effects which yield the black histograms plotted) the 
blue dashed line in the upper right panel. Although the agree- 
ment may appear reasonable, the difference is significant - 
such a fit overpredicts the fraction of high redshift objects 
at Eddington ratios < 0.1 and underpredicts the fraction at 
^0.3, giving a somewhat poor fit overall {x^ /v = 2.1, but with 
typical > 3a overpredictions for Eddington ratios < 0.1). 

Furthermore, without being modified to allow for a distribu- 
tion of Eddington ratios, such models are clearly inconsistent 
with the observations, as they would predict a single, constant 
Eddington ratio. However, models which fit the observed evo- 
lution in the quasar luminosity function with a non-static dis- 
tribution of accretion rates do recover the broadening of the 
Eddington ratio distribution at low redshift, so long as strong 
evolution in the distribution of accretion r ates for systems of 
a giv en black hole mass is not allowed JSteed & Weinberg 
giving a qualitatively similar picture of the evolution 
we model. Regardless of the evolution in accretion rates, an 
advantage of our modeling is that it provides a physically mo- 
tivated predicted distribution of accretion rates, as opposed to 
being forced to adopt the distribution of accretion rates from 
observational input (which can be, as demonstrated in the fig- 
ure, significantly biased by observational selection effects). 
The dotted histograms show the distribution if we apply our 
ADAF correction to the intrinsic distribution, and demon- 
strate that this does not significantly change the result. We 
note that our model for black hole accretion employs the Ed- 
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Fig. 2 1 . — Predicted distribution of Eddington ratios based on the luminos- 
ity function and the quasar evolution in our simulations, at three redshifts z = 
0.5 (top panels), z = 10 (middle), and z = 2.0 (bottom). The inferred distribu- 
tion of Eddington ratios, adopting a constant bolometric coirection from the 
observed (attenuated) luminosity in each of three bands is shown, i.e. assum- 
ing L = ULf'' (4400 A; left), L = 52L°^" (0.5-2 keV; middle), and L = 35 L°^^ 
(2-10 keV; right). For each waveband, results are shown for three reference 
luminosities. In B-band, Mg < -19 (red), Mg < -22 (blue), and Mb < -25 
(black). In soft X-rays, log(Lix[ergs"']) > 40 (red), 42 (blue), 44 (black). 
In hard X-rays, log(LHx[ergs"']) > 41 (red), 43 (blue), 45 (black). 



dington limit as a maximum accretion rate; if we remove this 
restriction, we find that the simulations spend some small but 
non-negligible time with / ^ 1 - 2, which is also consistent 
with the observations. 

Furthermore, we can make a prediction of this model which 
can be falsified, namely that the Eddington ratio distribution 
at luminosities well below the break in the luminosity func- 
tion should be broader and extend to lower values than the 
distribution at luminosities above the break luminosity. We 
quantify this in the lower panels of Figure |20| for the distri- 
bution at low redshifts z < 1 . Here we consider two bins in 
luminosity, L = lO^-^ - IO^'-^Lq and L = lO'^-^ - IO'^'^Lq, for 
redshifts where the break in the luminosity function is at ap- 
proximately L ^ 10" - IO'^Lq. Clearly, the distribution is 
broader and extends to lower Eddington ratios in the former 
luminosity interval, whereas in the latter it is strongly peaked 
about / ~ 0.2- 1, for both the complete distribution (black) 
and that with LB,ob,s > 10" (red). The distribution obtained 
applying the ADAF correction described above is shown as 
dotted histograms. Despite the fact that the Eddington ratio 
distribution at low luminosities will be strongly biased by se- 
lection effects, a reasonably complete sample should be able 
to test this prediction, at least qualitatively. 

We illustrate the effects of changing observed waveband, 
redshift, and luminosity thresholds on the observed Edding- 
ton ratio distribution in Figure]^ Here, we plot the pre- 
dicted distribution of Eddington ratios determined as in Fig- 
ure |20| from our fitted n(Lpeak) distribution at each redshift 
and the distribution of Eddington ratios as a function of in- 
stantaneous and peak luminosity in each of our simulations 
(specifically, these are drawn from the Monte Carlo real- 
izations of the quasar population described in § |8|i. We 
show the predictions at three redshifts z = 0.5 (top panels), 
z = 1.0 (middle), and z = 2.0 (bottom). For each redshift, 
results are shown in three wavebands, and with three refer- 
ence luminosities. In B-band, we require Mb < -19 (red), 
Mb < -22 (blue), and Mb < -25 (black). In soft X-rays, 
log(Lsx[ergs-']) > 40 (red), 42 (blue), 44 (black). In hard 
X-rays, log(LHx[ergs-i]) > 41 (red), 43 (blue), 45 (black). 
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The observationally inferred distribution of Eddington ratios 
at each redshift is loosely estimated by adopting a constant 
bolometric correction from the observed (attenuated) lumi- 
nosity in each of three bands shown, i.e. assuming L = \2U^^ 
(4400 A; left), L = 52L°\' (0.5-2 keV; middle), and L = 35Lg'| 
(2-10 keV; right). This follows common practice in many ob- 
servational estimates of the Eddington ratio distribution and 
allows for the effects of attenuation, but we caution that it can 
be misleading. 

If we instead use the lum inosity-dependent bolometric cor- 
rections of .Marconi et al.i C2004.) which we adopt through- 
out, even given that we are calculating from the observed 
(attenuated) luminosities, we do not see the large population 
of highly sub-Eddington (Eddington ratios < 10"^) quasars 
in soft and hard X-ray samples with low luminosity thresh- 
olds. This is because these are actually reasonably high- 
Eddington ratio quasars, but our bolometric corrections im- 
ply that a larger fraction of the bolometric luminosity is ra- 
diated in the X-ray at low bolometric luminosity, meaning 
that assuming a constant bolometric correction will under- 
estimate the Eddington ratios of high-bolometric luminosity 
sources. Regardless, the figure illustrates both the importance 
of different wavelengths (i.e. the ability to observe more low- 
Eddington ratio sources in X-ray as compared to optical sam- 
ples) and luminosity/magnitude limits on the inferred distri- 
bution of Eddington ratios. For example, even for relatively 
deep B-band quasar samples complete to Mb < -23 (i.e. com- 
plete to essentially all objects traditionally classified as having 
"quasar-like" luminosities), the expected observed Eddington 
ratio distribution at z ^ 0.5-2 is quite sharply peaked about 
~ 0.1 -0.3, in good agreement with recent observational re- 
sults (|Kollmeier et al. 2005). 

We do not compare to the z = Q distribution of black 
hole accretion rates, as this is dominated by objects at 
extr emely low Eddington ratios I ~ 10"^ - 10""^ (e.g., 
|HoI 12002 Marchesini et al. 200^ Uested l2005h . which 
are well below the range we model, and are not likely 
to be driven by merger activity (many of these objects 
are quiescent, low-luminosity Seyferts in normal spi- 
ral galaxy hosts); furthermore, many of these objects 
are not accreting at th e Bondi rate ( Fabian & Canizares. 
1988[ IBlandford & Beg elman 1999; Di Matteo et all 
200d 'Naravan et al.l '2000; Ouataert & Gruzinov 2000^ 
Di Matteo et al. 2001; Loewenstein et al. 2001; Bowe r et aP 
2003i) . clearly showing that our simulations must incorporate 
more sophisticated models for accretion in quiescent, low- 
luminosity states (when gravitational torques cannot provide 
a mechanism to drive large amounts of gas to the central 
regions of the galaxy) in order to describe such phases. 

However, it has been suggested that the rapid "blowout" 
phase and subsequent decay in accretion rates seen in our 
simulations, coupled with spectral modeling of radiatively 
inefficient accretion modes, can explain the apparently bi- 
modal distribution of low-redshift accretion rates (Cao & Xu 
l2005h . Moreover, present-day, relaxed ellipticals are observed 
to have mass accretion rates ^ 10"^ implying a long relaxation 
time at moderate and low accretion rates, qualitatively similar 
to that seen after the "blowout" in our modeling (Hopkins et 
al. 2005f). A pure exponential decay in accretion rate after the 
peak quasar phase would give m = M/MEdd ^ ^^Vi~tH /tg) at 
present, where tn is the Hubble time and fg is the quasar life- 
time of order e.g. the Salpeter time = 4 x 10^ yr, yielding an 
unreasonably low expected accretion rate m ^ lO"''*^. Even 



assuming an order of magnitude larger quasar lifetime, this 
gives rii ^ 10"'^, far below observed values, implying that re- 
gardless of the fueling mechanisms at low luminosities, the 
basic key point of our modeling must be true to some extent, 
namely that quasars spend long times relaxing at moderate to 
low Eddington ratios. 

6. THE MASS FUNCTION OF RELIC SUPERMASSIVE BLACK 
HOLES FROM QUASARS 

From the Mbh-c relation and other host galaxy-black 
hole scalings, estimates of bulge and spheroid veloc- 
ity dispersions have been used to determine the to- 
tal mass density (pbh) and mass distribution of lo- 
cal, primarily inactive supermassive black holes (e.g., 
'Salucci et al."'1999', 'Marconi & Salvati"2002'; 'Yu & Tremaind 
2002; Ferrarese 2002; Aller & Richstone 2002; Marconi et ay 
l2004t IShankar et alJ I2OO40 . These estimates, along 
with others based on X-ray background synthesis (e.g., 
Fabian & lwasawal ll999t IeMs et alJl2002h . have compared 
these quantities to those expected based on the mass distribu- 
tion of 'relic' black holes grown in quasars. It appears that 
most, and perhaps nearly all of the present-day black hole 
mass density was accumulated in bright quasar phases, and 
the Mbh - cr and Mbh - ^buige correlations yield estimates of 
the local mass function in good agreement with those from 
hard X-ray AGN luminosity functions (Marconi et al. 2004). 

However, this modeling is dependent on several assump- 
tions. Namely, the average radiative efficiency e^, Eddington 
ratio /, and average quasar lifetime fg are generally taken to be 
constants and either input into the model or constrained by de- 
manding agreement with the local mass function. In our sim- 
ulations, we find the quasar lifetime and Eddington ratio to be 
complex functions of both luminosity and host system proper- 
ties (as opposed to being constants). We also find that quasars 
spend a large fraction of their lives in obscured growth phases, 
suggesting some mass gain outside of the bright quasar phase. 
It is thus of interest to determine the relic black hole mass 
function expected from our model for quasar evolution. 

Using our estimate for the birthrate of quasars with a given 
peak luminosity at a particular redshift, «(Lpeak), obtained 
from the luminosity function in § 13.21 we can estimate the 
total number density of relic quasars accumulated by a partic- 
ular redshift that were born with a given Lpeak (per logarithmic 
interval in Lpeak) from 

«(W = / «(Wdf = / (34) 

By redshift z = 0, most of these quasars will be "dead," with 
only a small residual fraction having been activated in the re- 
cent past. 

Using our log-normal form for n(Lpeak), with normaUza- 
tion n, and dispersion ct* held constant and only the me- 
dian = exp(^ir) evolving with redshift, this integral can 
be evaluated numerically to give the space density of relic 

quasars «(Lpeak)- Finally, we use Mgjj(Lpeak), roughly the Ed- 
dington mass of th e giv en peak luminosity (but determined 
more precisely in § I2.4> to convert from dn(Lpeak)/dlogLpeak 
to dn(MBH)/dlogMBH. This formulation implicitly assumes 
that black holes do not undergo subsequent mergers after the 
initial quasar-producing event. However, this effect should 
be small (a factor < 2) as subsequent mergers would be dry 
(gas poor). We explicitly calculate the effects of dry merg- 
ers on the spheroid mass function (essentially a rescaling of 
the black hole mass function calculated here) in lHopkins et alJ 
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and show that this is a small effect (significantly less 
than the uncertainties owing to our fit to the quasar luminos- 
ity function) even assuming the maximum dry merger rates of 
e.g.|van Dokkum (2005). 

This mass function can then be integrated over c/Mbh to 
give the total present-day black hole mass density, pbh- Ne- 
glecting temporarily the small corrections to Mgjj(Lpeak) from 
§ 12.41 we expect 

^peak is 



^BH ~ ^Edd(ipeak) = 



(35) 



where ts/eyC^ w 2.95 x 10 ^Mq /Lq, so therefore, 
ts f 

PBH= 2 / ^peakn(i'peak)dlogLpeak- (36) 

This can be combined with the integral over redshift for 
n(Lpeak), giving, at each z, a pure Gaussian integral over 
log(Lpeak), in the form 



Pbh- 



i5_^gi(<T.lnlO)- / 



(l+z)H(z) 



(37) 



where H(z) = H{z)/Hq and ry is the fractional lookback time 
at some upper limit. We must modify this integral above 
z ~ 2 to account for the decreasing space density of bright 
quasars, applying eith er ou r density or peak-luminosity evo- 
lution turnover from § 13.21 but quasars at these high redshifts 
contribute only a small fraction to the present-day density. 
Thus, in this formulation, the evolution of the total supermas- 
sive black hole mass density, i.e. pmi{z) / p^wiz = 0), is given 
approximately by the dimensionless integral above, and de- 
pends only on how evolves, essentially the rate at which 
the break in the quasar luminosity function shifts. Although 
this is not strictly true if we include corrections to Mgj^(Lpeak) 
based on Lpeak, the difference is small and this behavior is es- 
sentially preserved. Note that the total supermassive black 
hole mass density is independent of corrections from sub- 
sequent dry mergers, which (being gas poor) conserve total 
black hol e m ass. 

Figure |22] shows our prediction for the mass distribution 
of supermassive black holes, as well as the total density 
Pbh and its evolution with redshift. We find a total relic 
black hole mass density of pbh = 2.9'!i\ l ^ lO^MoMpc"-', in 
agreement with the observational estimate o f pbh = 2.9 ± 
0.5 hlj X 10^MoMpc-^ by Yu & Tremain^ (1200211) (hoj = 
//o/70kms~' Mpc~'; their result is converted from h = 0.65), 
an d within la of t h e value Pbh = 4.6t\ lhlj x 10^ Mpc"^, 
of 'Marconi et al. ' ("^OP), based on the observations of 
^ai-zke et al. ( 1994), Kochanek et al. (2001 ), Nakamura et al. 
( 120031) . iBernardi et alJ (I2003t) . and Sheth et al. (2003). The 
fractional evolution of pbh with redshift is quite well con- 
strained, and we find, as with previous estimates, that most of 
the present-day black hole mass density accumulates at mod- 
erate to low redshifts z ~ 0.5-2.5. The la errors are shown 
as dotted lines in the figure, and are close to our best-fit es- 
timate, as we have demonstrated that this quantity depends 
only on ki, the rate of evolution of the break in the lumi- 
nosity function with redshift, which is fairly well-constrained 
by observations (from our fitting to the luminosity functions, 
A:/, = 5.61 ± 0.28). The difference in pbh if we include or ne- 
glect the small corrections to is negligible compared to 
our errors 5%). 



Our estimate for the relic black hole mass distribution (thick 
black line) also agrees well with observational estimates, 
with all observations within the range allowed by the la er- 
rors of our fitting to the luminosity fun ction (dotted lines). 
The observations shown are again from Marconi et al. (2004|, 
based on the combination of observations by Marzke et ay 
(1994), Kochanek et al. (2001), Nakamura et al. (2003), 
iBernardi et aL ti2003i) . and iSheth et al.i(>2003t) . The high mass 
end of the black hole mass function Mbh > 10^ is rel- 
atively sensitive to whether or not we apply the Mgjj(Lpeak) 

corrections of § 12.41 instead of taking Mgjj = MEdd(i'peak) (thin 
line), as well as to our fitting procedure. However, the agree- 
ment is still good, and this is also where the observational es- 
timates of the mass distribution are most uncertain, as they are 
generally extrapolated to these masses, and are sensitive to the 
assumed intrinsic dispersions in the Mbh-t and MBH-^buige 
relations ( Yu & Tremaine 2002). 

If, instead, we adopt a light-bulb, constant Eddington 
ratio, or exponential light curve model for quasar evolu- 

tion, we would have Mgjj oc Lpeak, and thus the prediction 
would be similar to the thin black line shown, a some- 
what worse fit at high black hole masses. However, in 
these models this can be remedied by adjusting the typi- 
cal Eddington ratios, quasar lifetimes, or radiative efficien- 
cies. We do not show the range of predictions of these 
models for the relic supermassive black hole mass func- 
tion, as they have been examined in detail previously (e.g. 
Salucci et al. 1999; Marconi & Salvati 2002; Yu & TremainJ 
.2002j .Ferrarese 2002; Alle r & Richstone.2002; Maixoni et aL 
l2004t lShankar et al. 2004). These works demonstrate that the 
observed quasar luminosity functions are consistent with the 
relic supermassive black hole mass function, given typical ra- 
diative efficiencies e,- 0.1 and Eddington ratios ^ 0.5- 1.0, 
and that most of the mass of black holes is accumulated in 
bright, observed phases, or else the required radiative effi- 
ciency would violate theoretical limits. 

That our model of quasar lifetimes and obscuration repro- 
duces the observed z = supermassive black hole mass func- 
tion explicitly demonstrates that we are consistent with these 
constraints. By choice, the radiative efficiency in our simula- 
tions is = 0.1, and accretion rates are not allowed to exceed 
Eddington. As noted in §|3 most of the black hole mass is ac- 
cumulated and radiant energy released in the final, "blowout" 
phase of quasar evolution, and here our black hole mass func- 
tion and cumulative black hole mass density demonstrate that 
our modeling is consistent with integrated energy and mass 
arguments such as that of Soltan ( 1982), despite the fact that 
quasars spend more time in obscured phases than they do 
in bright optical quasar phases. In fact, comparison of our 
predicted total black hole mass density with estimates from 
the z = black hole mass distribution allows some latitude 
for significant mass gain in radiatively inefficient growth or 
black holes in small, disky spheroids, although we emphasize 
that this is mainly because the uncertainty in our prediction is 
large, it is not inherent or necessary in our modeling. 

The anti-hierarchical nature of black hole formation, where 
less massive black holes are formed at lower redshift, is re- 
flected in our modeling by the shift of the break in the quasar 
luminosity function to lower values with decreasing redshift. 
This can be seen in Figure |22] where the black hole mass 
distributions are shown at redshifts z = 1.5, 3.0 and 5.0, as- 
suming either pure peak luminosity evolution or pure den- 
sity evolution for z > 2 (dot-dashed and dashed, respectively). 
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Fig. 22. — Right: Total predicted quasar relic black hole mass density and evolution of the fractional black hole mass density with redshift. Dotted lines show 
the difference resulting from Icr deviation in fitted H(Lpeak) from the luminosity function. Left: Predicte d present z = relic mass function (thick black line), for 
comparison with the Ic range (yellow) of the infeired supermassive black hole mass function from .Marconi et all 120040 . Also shown are the results given Icr 
eiTors in the fitted n(Lpeak) disfiibution (dotted lines), or ignoring the small connections to M^^iLp^..^^^) from § l2.4l (thin black line). Dot-dashed lines show the 
predicted mass function at z = 1.5, 3.0, 5.0 (blue, green, and red, respectively). The extensions to z > 2 includes the turnover (pure peak luminosity evolution 
fonii) in the quasar space density above z = 2 from high-redshift luminosity functions described in § 13.21 except for the dashed green and red lines which use the 
pure density evolution form. 
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Fig. 23. — Fractional number density n{M, z)/n{M, z = 0) of black holes of a given mass as a function of redshift, for several different black hole masses 
as shown. Fo r z > 2 this includes the turnover (pure density evolution form) in the quasar space density above z = 2 from high-redshift luminosity functions 
described in § 13.21 Left panel shows the results using our full model of quasar lifetimes, right panel assuming a "light-bulb" or exponential (constant Eddington 
ratio) light curve model. The yellow dot-dash (IO'Mq) and red triple-dot-dash (IO'^Mq) curves are nearly identical in the right panel. 



While the choice for the turnover in the z > 2 quasar den- 
sity matters little for the z < 2 black hole mass functions, 
the Iow-Mbh distribution at high redshift (where observations 
do not constrain n(Lpeak) well) is quite different between the 
two models. Figure |23l plots the fractional number density 
of black holes of a given mass as a function of redshift, i.e. 
n(M, z)/n(M,z = 0), where n(M) = dn/dlog(M) is just the 
number density at mass M. This figure demonstrates that 
higher-mass black holes originated over a larger range of red- 
shifts, and that they mostly formed at higher redshift, com- 



pared to lower-mass black holes. 

The right panel of Figure |23l compares our prediction 
to that of a light-bulb or exponential light curve model 
for quasar lifetimes. In these models, the anti-hierarchical 
nature of black hole assembly is dramatically suppressed. 
At the high-mass end, there is no measurable difference 
in the distribution of formation redshifts (i.e. the Mbh = 
IO'^Mq and Mbh = 10"' Mq curves are indistinguishable), 
and there is little change in the formation times at Mbh = 
IO^Mq. The shift in formation redshift at lower masses. 
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although significant, is smaller than that predicted in our 
model. If spheroids and black holes are produced to- 
gether, as in our picture, these models of the quasar life- 
time would imply that spheroids of masses Myir ^ 10"- 
10'^ Mq all formed over nearly identical ranges of redshifts, 
which is inconsistent with many observations indicating anti- 
hierar chical groy yth of the red, elliptical galaxy population 
(e.g. iTreu et al.' '2001', 'van Dokkum et al.' '2001', Treu et al.' 
2002t van Dokkum & Stanford 2003; Gebhardt et al. 2003; 
Rusin et al. 2003; van de Ven et al. 2003; Wuvts et al. 2004t 
Treu et al. 2005; Holden et al. 2005; van der Wei et al. 20051 
di Serego Alighieri et al. 2005; Nelan et al. 2005). Implica- 
tions of our model for the red galaxy sequence are considered 
in Hopkins et al. (2005e), where we show that this weaker 
anti-hierarchical black hole (and correspondingly, spheroid) 
evolution is inconsistent with observed luminosity functions, 
color-magnitude relations, and mass-to-light ratios of ellipti- 
cal galaxies. 

Our modeling reproduces the observed total density and 
mass distribution of supermassive black holes at z = with 
black holes accreting at the canonical efficiency e, =0.1 ex- 
pected for efficient accretion through a Shakura & Sunvae_3 
111973) disk. Presumably, a large change in e, would give a 
significantly different relation between peak luminosity and 

black hole mass (for the same Lpeak, -^bh oc 1 / f r), and thus 
if the quasar lifetime remained similar as a function of peak 
luminosity, this would translate to a shift in the black hole 
mass function. The long obscured stage in black hole evolu- 
tion does not generate problems in reproducing the black hole 
mass density, and the final phases of growth are still in bright 
optical quasar stages. However, a large Compton-thick popu- 
lation of black holes at all luminosities (or even at some range 
of luminosities at or above the break in the luminosity func- 
tion) (e.g., Gilli et al. 2001; Ueda et al. 2003), or a large pop- 
ulation accreting in a radiatively inefficient ADAF-type so- 
lution, as invoked to explain discrepancies in the X-ray back- 
ground produced by synthesis models (Di Matteo et al. 1999), 
would result in a significant over-prediction of the present- 
d ay su permassive black hole density. As we demonstrate in 
§ 17.21 invoking such populations is unnecessary, as our pic- 
ture for quasar lifetimes and evolutionary obscuration self- 
consistently reproduces the observed X-ray background. 

Finally, we note that we reproduce the z = distri- 
bution of black hole masses inferred from the di stribu- 
tion of spheroid velocity dispersions (Sheth et al. 2003) and 
luminosity functions (Marzke et al. 1994; Kochanek et aTl 
1^01; Nakamuraet al. 20()3|), based on the observed Mbh- 
cr relation and fundamental plane for galax y properties 
(e.g.,[Bernardietal. 2003; Gebhardt etan i2003h . Therefore, 
since our modeling also reproduces the observed Mbh - <J 
(|pi Matteo et al. 2005; Robertson et al. 2005b) and funda- 
mental plane (Robertson et al., in preparation) relations, we 
implicitly reproduce the z = distribution of spheroid velocity 
dispersions and spheroid luminosity functions, given our ba- 
sic assumption that the mergers that produce these spheroids 
also give rise to luminous quasar activity. 

7. THE COSMIC X-RAY BACKGROUND 

7.1. The Integrated Spectra of Individual Quasars 

Unresolved extragalactic sources, specifically obscured 
AGN, have been invo ked to explain the cosmic X-ray back- 
ground (e.g, ISetti & Woltier 1989) . This picture has been 
confirmed as deep surveys with Chandra and XMM-Newton 



have resolved most or all of the X-ray background into 
discrete sources, p r imarily obscured and unobscured AGN 
(iBrandt et al."'2001', 'Hasinger et al."2001Ji iRosati et al.ll2002t 
Giacconi et al. 2002; Baldi et al. 2002). The X-ray back- 
ground, however, has a harder X-ray spectrum than typical 
quasars, with a photon index F ~ 1.4 in the 1 - lOkeV range 
(jMarshall et al. 1980). Therefore, obscured AGN are impor- 
tant in producing this shape, as absorption in the ultraviolent 
and soft X-rays hardens the observed spectrum. Indeed, pop- 
ulation synthesis models based on observed quasar luminosity 
functions and involving large numbers of obscured AGN have 
been successful at matchi ng both the X-ray background in- 
tensity and spectral shape (Ma dau et all 19941 IComastri et al.l 
1995; Gilli et al. 1999, 2001). However, these models make 
arbitrary assumptions about the ratio of obscured to unob- 
scured sources and its evolution with redshift, choosing these 
quantities to reproduce the X-ray background. Furthermore, 
as X-ray surveys have been extended to higher redshifts, it 
has become clear that both the observed redshift distribution 
of X-ray sources and the ratio of obscured to unobscured 
sources is inconsistent with that required by these models 
(|Hasinger 2002; Bargeretal. 2003). Even synthesis models 
based on higher-redshift X-ray surveys and using observation- 
ally derived rati os of obscured to unobscured sources (e.g., 
Ueda et al. 2003) have invoked ad hoc assumptions about ad- 
ditional populations of obscured sources to reproduce the X- 
ray background shape and intensity. 

We can test our model by examining whether the quasar 
luminosity function, relic AGN mass distribution, and X- 
ray background can be simultaneously reproduced in a self- 
consistent manner Because our formulation describes the 
birthrate of quasars with a peak luminosity Lpeak, it is most 
useful to consider the integrated energy spectrum of such a 
quasar over its lifetime, 

vE, = jdtvL,{t) = j vUiDL^^^^^dlogL, (38) 

where fi,{L) is the bolometric correction {Li, = f^L). As 
an example. Figure |23 shows the integrated intrinsic spec- 
tra (thick solid lines) from the simulations Al, A2, A3, A4, 
and A5, described in § 12.11 The final black hole masses for 
these simulations are M^^^ = 7x 10*, 3x 10'', 3 x 10^, 7x 
10**, and 2 x \(f Mq, respectively. The integrated spectral 
shape in the X-ray, in particular, is ultimately determined 
by the observationally motivated bolometric corrections of 
IMar coni et al. ( 2004), with a reflection component in the 
X-ray determined following Magdziarz & Zdziarski ( 1995), 
and, in the case of the observed spectrum, the distribution of 
column densities calculated from the simulations. Using our 
fits to the lifetime df/dlogLasa function of instantaneous and 
peak luminosities, we can calculate the expected vEy from 
the integral above. These integrated spectra are shown as the 
dot-dashed lines in the figure, and agree well with the actual 
integrated spectra of the simulations, demonstrating the self- 
consistency of our model and applicability of our fitted life- 
times. 

This can be compared to idealized models for the quasar 
lifetime, where we allow the quasar to radiate just at its peak 
luminosity Lpeak ~ LEdd(A^BH) f^^" some fixed lifetime fg. We 
determine fg by demanding that the total energetics be correct, 

Lpeakfg = ^rM^^c^. The predicted integrated energy spectra 
are shown as the dashed lines, and under-predict the soft and 
hard X-ray energy output by a factor ^ 1.5-2. This is be- 
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Fig. 24. — Left: Integrated intrinsic spectra (tliiclc solid lines) from simulations Al, A2, A3, A4, and A5 (black, blue, green, yellow, red, respectively), with 
virial velocities Vyji = 80, 113, 160, 226, and 320 kms"' . The predicted integrated spectra from our model for quasar lifetimes are shown as dot-dashed lines, and 
the prediction of a "light bulb" model, where the same total energy is radiated at L = Lpe^ki as dashed lines. Integrated observed spectra are shown as thin solid 
lines. Right: Integrated observed X-ray spectrum from the A3 s imul ation (thick black line), compared with the integrated intrinsic spectrum, reddened by various 
column density distributions: our fitted Nn distributions from ii l2.3l (thick black dashed Une), constant (luminosity-independent) lognormal A'h distribution with 
Nh = 10-^ cm"- and aNf, = 0.4, 0.7, 1.0 (blue, green, and red dashed lines, respectively), and constant A'h = 10^', 10^'^, 10^^, 10^^'^, and 10-^ cm"- (thin 
dot-dashed lines). 



cause higher-luminosity quasars tend to have a larger fraction 
of their e nergy radiated in the UV-optical rather than the X- 
ray (e.g.. •Wilkes et al.""1994'. 'Gree n et alJlT99l iVignah et alJ 
|2()03; Strateva et al. 2005), reflected in our bolometric cor- 
rections. Thus, assuming that the quasar spends all its time 
at ipeak does not account for extended times at lower lumi- 
nosity, where the ratio of X-ray to total luminosity is higher, 
which would generate an integrated spectrum with a larger 
fraction of its energy in the X-ray. Assuming that the quasar 
undergoes pure Eddington-limited growth to its peak luminos- 
ity produces an almost identical integrated spectrum to this 
light-bulb model, as it is similarly dominated by L ~ i-peak- 

Of course, the intrinsic integrated energy spectrum of the 
simulations is not what determines the X-ray background, but 
rather the integrated observed spectrum is the critical quantity. 
This is shown as the thin lines in the left panel of Figure l24l 
and in detail for our fiducial A3 simulation in the right panel 
of the figure (thick solid line). Along a given sightline, the 
observed integrated spectrum will be 

where t^, is the optical depth at a given frequency. We can 
integrate over solid angle and obtain 

^£..obs = / ^/.(.--)L^^^^^dlogL, (40) 
J dlogL 

where (e~^'^) is the averaged e"''" over the column den- 
sity distribution P{NH\L,Lj,eai_). Using our fits to the col- 
umn density distribution and quasar lifetimes and calcu- 
lating i/E ^ obfi as above, we reproduce the integrated ob- 
served spectrum quite well (black dashed line). For com- 
parison, we show that it is not a good approximation to red- 
den the spectrum with a constant A^h, giving the results for 
A^H = lO^i, 10-' ■^ 10^2, 1022•^ and lO^^ cm'^ (thin dot-dashed 



lines). Even allowing for a distribution of A^h values, the re- 
sulting spectrum is a poor match to the observed one if that 
distribution is taken to be static (i.e. luminosity-independent, 
as in traditional torus models, for example). We show the 
results of reddening the intrinsic spectrum by such a (Gaus- 
sian) distribution, varying the dispersion (Tnh = 0.4,0.7, 1.0 
(blue, green, and red dashed lines, respectively), for a median 
column density Nh = 10^^ cm~^, the median column density 
expected around Lpeak in this simulation. Therefore, the lumi- 
nosity and host system property dependence of both quasar 
lifetimes and the column density distribution must be ac- 
counted for in attempting to properly predict the X-ray back- 
ground spectrum from observations of the quasar luminosity 
function. Finally, note that the hard cutoff in the observed UV 
spectra at 912A owes to our calculated cross-sections being 
incomplete in the extreme UV. Properly modeling the escape 
fraction and observed emission at these frequencies, while not 
important for the X-ray background, is critical to calculat- 
ing the contribution of quasars to reionization, and requires 
a more detailed modeling of scattering and absorption, espe- 
cially in the bright optical quasar phase. 

7.2. The Integrated X-Ray Background 

Given the volume emissivity j,^(z) (per unit comoving vol- 
ume) of some isotropic process at a given frequency at redshift 
z, the resulting background specific intensity at frequency V() 
atz = is (Peacock 1999) 

If we were to consider the emissivity per unit physical vol- 
ume, there would be an extra factor of ( 1 -I- z)~^ in the integral 
above. In § 17.11 we determined the integrated observed en- 
ergy £'i/.obs(i'peak) produccd by a quasar with peak luminosity 
i'peak- We have also inferred n(Lpeak)(z) in § 13.21 the rate at 
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which quasars of peak luminosity Lpeak are created per unit 
comoving volume per unit cosmological time. Therefore, the 
comoving volume emissivity is just 

peak 5 

(42) 

or, expanding E^^obs, 

jAz) = y^/logLpeak J dlogL 

X/.(.--)L^^5^^«(Lpeak). (43) 

dlogL 

If the column density distribution were independent of Lpeak. 
as is assumed in even luminosity-dependent torus models or 
observationally deter mined A^h functions used for X-ray back- 
ground synthesis (e.g.. lUeda et al..2003.) . then we could com- 
bine terms in Lpeak and integrate over them. This simplifica- 
tion, along with the definition of the luminosity function in 
terms of Lpeak, gives the more traditional formula for the X- 
ray background in terms of only the observed column density 
distribution and luminosity function, 

jAz)= dlogL— —-L,{e-^^). (44) 
J dlogL 

However, as we showed in § 13.31 and § 17.11 neglecting the 
dependence on Lpeak is not a good approximation at all lu- 
minosities and gives an inaccurate estimate of the integrated 
quasar spectrum; therefore, "purely observation-based" syn- 
thesis models of the X-ray background will be inaccurate in 
a similar manner to synthesis models with an inappropriate 
model for the quasar lifetime. Essentially, this "averages out" 
the varying distribution of column densities with Lpeak, which 
changes the shape of the spectrum in a non-linear manner, es- 
pecially when integrated over varying bolometric corrections 
as shown above. 

Figure |25l (upper panel) shows the predicted X-ray back- 
ground spectrum from our full modeling of quasar lifetimes 
and obscuration (solid lines). We use our analytical fits to the 
quasar lifetime and column density distributions as in § 17.11 
above, as Figure |24] demonstrates that they accurately repro- 
duce the actual integrated quasar X-ray spectra of the simu- 
lations, and the analytical forms are integrated over all lumi- 
nosities and redshifts. The dotted lines show the deviation 
resulting from shifting the parameters describing our fitted 
n(Lpeak) distribution by 1 cr in either direction, although degen- 
eracies in the parameters suggest that the actual uncertainty in 
the background prediction is smaller The dashed line shows 
the predicted X-ray background if we ignore the broadening 
of the A^H distribution across simulations (ctaj^ = 1 .2) and in- 
stead consider only the dispersion of an individual simulation 
at a given luminosity ((Ta(„ = 0.4). 

These can be compared to the observations of Gruber et al. 
^39% (red curve, for E > 3keV) and Barcons et al. ( 2000.) 
(cyan curve, for E < lOke V). We increase the normalization 
of the lGrub er et al. spectrum to match that of the best 

estimate from Barcons et al. ( 2000) over the range of overlap, 
determined from combined ASCA, BeppoSAX, and ROSAT 
data to be 10.0:l:o 9keV cm"^ s"' sr"' keV"' at 1 keV. The uncer- 
tainty in the normalization between the two samples, 20%, 
is shown as the shaded yellow range (alternatively, this repre- 
sents the ~ 2cr errors in the /?(9&4r normalization). 

In the middle panel of the figure, we calculate the predicted 
X-ray background using our full model of the quasar life- 
time, but with different models for quasar obscuration. The 
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Fig. 25. — Predicted integrated X-ray background spectrum (solid black 
line) from our model of quasar lifetimes and attenuation, with the peak lu- 
minosity distribution n(Lpeak) determined from the luminosity function. Blue 
and red thick line s show the observed spectrum from Barcons et al. 1 2()0(|) 
and Gru ber et al J tT99 9 ). respectively. The shaded yellow area illustrates the 
uncertainty in normalization between both samples (alternatively, la errors 
in the Barcons et al. 2000 normalization). The predictions given Icr devia- 
tions in the fitted «(Lpeak) distribution (dotted lines) and given the n(Lpeaii) 
distribution determined from hard X-ray data only (dashed line) are shown in 
the upper panel. Middle panel shows the prediction using our modeling of 
quasar lifetimes but different models of obscuration, lower panel the predic- 
tion with a "light-bulb" or exponential (constant Eddington ratio) model and 
different obscuration models. 



solid black line shows the prediction using our full model 
of quasar obscuration, and is identical to the solid black line 
in the upper panel. The observations are likewise shown in 
an identical manner to the upper panel. The dashed black 
line is the prediction adopting the standard torus model for 
quasar obscuration, and the dotted line adopts the receding 
(luminosity-dependent) torus model. These models produce 
the same overall ^ 30keV normalization, as this is relatively 
unaffected by obscuration, but they predict a slightly (^ 20%) 
higher background at low energies, giving a slightly softer 
spectrum. This may appear counterintuitive, given that in 
Figure these models tend to overpredict the number of 
high-column density sources, but this is because these models 
predict a strongly bimodal column density distribution, with 
unobscured sightlines encountering negligible column densi- 
ties. These unobscured sightlines dominate the soft X-ray in- 
tegrated spectrum, where the large column densities through 
the torus attenuate the quasar spectrum heavily. However, 
this net offset in the predicted background spectrum is gen- 
erally within the range of the systematic theoretical and ob- 
servational uncertainties, and can further be alleviated by tun- 
ing the parameters of the torus model to fit the X-ray back- 
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ground spectrum (e.g. Treister & Urry 2005, although their 
fits require a larger fraction of Compton-thick A^h ~ 1 0^^ cm"^ 
sources than shown for even the receding torus model in 
Figure \T2\. The feature at < 5keV in the standard torus 
model prediction is a consequence of assuming that "unob- 
scured" lines of sight encounter negligible column density, 
and does not appear if such sightlines encounter moderate 
(~ lO^'cm"^) columns. 

The lower panel of the figure shows the predicted X-ray 
background spectrum if we instead consider a light-bulb or 
exponential light curve (fixed Eddington ratio) model for the 
quasar lifetime, again with various descriptions of quasar ob- 
scuration. In such models, the predicted X-ray background 
spectrum is independent of the quasar lifetime or character- 
istic Eddington ratio assumed (see Equation I44>. However, 
as shown in Figure|23 these models do imply a different inte- 
grated spectrum for quasars; i.e. different effective bolometric 
corrections for predicting the X-ray background. In particular, 
in this model, the observed quasar spectrum at a given lumi- 
nosity (averaged over the quasar population at that luminos- 
ity) is the same as the "effective" quasar spectrum one would 
use to calculate the total contribution to the X-ray background 
from quasars of the corresponding observed or peak luminos- 
ity, whereas this is not true in our model of quasar lifetimes. 
The observations are shown in the same manner as the pre- 
ceding panels. The black solid line shows the prediction with 
this simplified model for the quasar lifetime, but still adopting 
our full model for obscuration as a function of instantaneous 
and peak luminosity, the dashed line assumes instead a stan- 
dard torus model for obscuration, and the dotted line assumes 
a receding torus for the obscuration. The variations among 
different obscuration models are relatively small at most en- 
ergies, and similar to those discussed above adopting our full 
model of quasar lifetimes. 

In all three cases, however, this model for the quasar life- 
time significantly under-predicts the X-ray background, par- 
ticularly at the ^ 30keV peak . This shortfall is well-known, 
and earlier attempts (e.g., Madau et alJll994l: IComastri et al. 
p95; GilH et al. 1999, 2001 : Pomo mo et alJ2000HUeda et al. 
l2003ft have generally had to invoke additional assumptions 
about large obscured populations or a strong increase in the 
obscured fraction with redshift, neither of wh ich is consistent 
with obser vations (e.g., Hasinger 2002; Barser et al."20Q3l 
lUeda et alJ 2003; Szoko lv etalJ l2004: Bargeretal. 200§. 
The difference between the predictions of various quasar life- 
time models is, as explained above, attributable to the differ- 
ence between the integrated quasar spectrum produced in our 
full model of the quasar lifetime (in which quasars spend long 
periods of time at low luminosities, with harder X-ray spec- 
tra), and the integrated spectrum in these simplified quasar 
lifetime models, which is proportional to the instantaneous 
quasar spectrum, and therefore underpredicts the hard X-ray 
portion of the spectrum by as much as ^ 50%. 

Our prediction of the X-ray background agrees well with 
the observed spectrum over the range ^ 1 - 100 keV. (At en- 
ergies above lOOkeV it is likely that processes we have not 
included, such as those involving magnetic fields, contribute 
significantly to the background.) Unlike previous synthesis 
models for the X-ray background, we are able to do so without 
invoking assumptions about large Compton thick populations 
or larger obscured populations at different redshifts. In part, 
this is because our modeling allows us to predict, based on 
n(i'peak) and our column density formulation, the population 
of Compton thick sources (see Figure fT2l . However, as we 



have demonstrated, it is primarily because the deficit in pre- 
vious synthesis models can be attributed to their inability to 
properly account for the dependence of quasar lifetimes and 
attenuation on both the instantaneous quasar luminosity and 
the host system properties (peak luminosity). Our picture, on 
the other hand, yields an estimate for the X-ray background 
spectrum that is simultaneously consistent with the observed 
supermassive black hole mass distribution and total density, 
as well as the "luminosity -dependent density ev olution" ob- 
served in X-ray samples (Ho pkins et alJl2005fh . The back- 
ground is primarily built up from z ~ 2.5 to z ~ 0.5, as is 
evident from the evolution of the black hole mass density in 
Figure|22l although a harder spectrum at low luminosities will 
weight this slightly towards lower redshifts (where more low- 
luminosity quasars are forming). Compton thick and relaxing, 
low-luminosity sources are accounted for, not as large, inde- 
pendent populations, but as evolutionary phenomena continu- 
ously connected to the "normal" quasar population. 

8. DISCUSSION 
8.1. General Implications of our Model 

Our modeling suggests two important paradigm shifts in 
interpreting quasar populations and evolution: 

(1) First, as proposed in Hopkins et al. ( 200^3), a proper ac- 
counting of the luminosity dependence of quasar lifetimes (as 
opposed to models in which quasars grow in a pure exponen- 
tial fashion or turn on and off as "light bulbs") implies a novel 
interpretation of the luminosity function. The steep bright 
end (luminosities above the "break" in the luminosity func- 
tion) consists of quasars radiating near their Eddington limits 
and is directly related to the distribution of intrinsic peak lu- 
minosities (or final black hole masses) as has been assumed 
previously. However, the shallow, faint end of the luminosity 
function describes black holes either growing in early stages 
of activity or in extended, quiescent states going into or com- 
ing out of a peak bright quasar phase, with Eddington ratios 
generally between / ~ 0.01 and 1. The "break" luminosity in 
the luminosity function corresponds directly to the peak in the 
birthrate of quasars as a function of peak luminosity n(Lpeak)- 

This interpretation resolves inconsistencies in a num- 
ber of previous theoretical studies. For example, semi- 
analy tical models of the quasar luminosity fun ctions 
(e,g. , 'Kauffmann & Haehneltll2000t iHaiman & Menoul»2000.: 
IWvithe & Loeb 2003) assume, based on simplified models for 
the quasar lifetime, that quasars at the faint end of the lumi- 
nosity function correspond to low final-mass black holes (low 
i'peak ~ L), presumably in small halos. Consequently, these 
models overpredict the number of active low-mass black holes 
(as estimated from radio sourc e counts), especially at high 
redshift, by orders of magnitude ('Haima n. Ouataert. & Boweg 
2004), and overpredict the number of low-mass spheroids and 
red galaxies observed dHopkins et al. 2005e'). 

Moreover, both observations iMcLure & Dunlop'20041) and 
comparison of the present-day black hole mass function with 
radio and X-ray luminosity functions (e.g. Merloni 2004) sug- 
gest anti-hierarchical evolution for the growth of supermas- 
sive black holes, where the most massive black holes were 
produced mainly at high (z > 2) redshift, and low-mass black 
holes mostly formed later, which does not follow from ideal- 
ized descriptions of quasar lifetimes and the luminosity func- 
tion (for a review, see e.g. Combes 2005). 

A one-to-one correspondence between observed luminos- 
ity and black hole mass does produce anti-hierarchical be- 
havior in some sense at the high-mass end, because the most 



Quasar Origins & Evolution 



39 



massive black holes are formed at z ~ 2 - 3 during the peak 
of bright quasar activity and the quasar luminosity function 
evolves to lower luminosities at lower redshifts (as is also 
the case for our model because the bright end of the lumi- 
nosity function is dominated by sources near their peak lu- 
minosities). However, at black hole masses equal to or be- 
low 10^ Mq (i.e. galaxies of stellar mass < IO^'Mq), the 
evolution in the quasar luminosity function implies a roughly 
constant production of black holes with these masses at all 
redshifts, which is inconsistent with observations of galaxy 
spheroids indicating that typical ages increase with mass, 
ruling out a large population of low-mass spheroids with 
ages equal to or older than thos e of hig h-mass spheroids 
(e.g.,_ p'reu et al l 1200 U Iv an Dokku m et 311120011 iTreu et al.l 



2002; van Dokkum & Stanford 2003; Gebhardt et al. 2003; 
Rusin et al . 2003; van de Ven et al. 2003; Wuvts et al. 2004; 
Treu et al.l!2 005; Holden et al. 2005; van der Wei et al. 2005t 



di Sereeo Alighieri et al..2005.: JMelan et al..2005.) . As demon- 



strated in Figure |23| such a model does not produce anti 
hierarchical growth or any age gradients within the high- 
mass spheroid population, also inconsistent with observa- 
tions. Even giv en observed "luminosity-dependent density 
evolution" (e.g. 'Pa ge et al.1 fT997t iMivaii et alJ I200a 
La Franca e t al. 200 21 ICowie et all l2003t lUeda et alJ 
etal. 2003; Hunt et al." '2004: 'Ciras uolo et al.1 



2001 



2003 



2005 



Fiore i 

Hasinger, Mivaii, & Schmidt 2005), implying that the densi- 
ties of lower redshift quasars peak at lower redshift, the in- 
ferred anti-hierarchical evolution if observed luminosity di- 
rectly corresponds to black hole mass (i.e. as in "light-bulb" 
or "fixed Eddington ratio" models) is not strong enough to ac- 
count for observed anti-hierarchical growth of the correspond- 
ing galaxy spheroids (Hopkins et al. 2005e). 

Furthermore, in these earlier models, a "break" in th e lumi- 
nosity function is not necessarily reproduced ( Wvithe & Loebl 
120031) . and the faint-end slope has no direct physical moti- 
vation. The break may be caused by feedback mechanisms 
which set a characteristic turnover in both the galaxy mass 
function and quasar luminosity function (e.g., Scannapieco & 
Oh 2004; Dekel & Birnboim 2004), as in our modeling. The 
n(Lpeak) distributions in our model and "light bulb" or "fixed 
Eddington ratio" models are comparable at and above the 
break in the quasar luminosity function, and therefore make 
similar predictions for some observations at these luminosi- 
ties. However, the faint-end slope has a different physical mo- 
tivation in our model. Unlike the bright-end slope, which is 
determined directly by the active final black hole mass func- 
tion or peak luminosity distribution (in essentially all models 
of the quasar lifetime), the faint-end slope in our modeling 
is a consequence of the quasar lifetime as a function of lu- 
minosity, and is a prediction of our simulations and model- 
ing almost independent of the underlying faint-end slope of 
the active black hole mass function or peak luminosity dis- 
tribution. In Ho pkins et al. (2005f ) we examine this in more 
detail, and demonstrate that it predicts well the evolution in 
the faint-end quasar luminosity function slope with redshift 
and the observed "luminosity-dependent density evolution" 
in many samples (Page et al. 1997; Mivaii et al. 2000, 2001'; 
l Ea Franca et al. 2002: Cowie et al. 2003; Ueda et al. 2003.; 
Fiore et all l2003riHunt et all 120041; ICirasuolo et alj bool 
iHasinger. Mivaii. & Schmidi2005h 

Other observational evidence for our picture exists; for ex- 
ample in the observed distribution of Eddington ratios (see 
§13), t he di stribution of low -redshift, active black hole masses 
(see § 14. 3> . and the turnover in the expected distribution of 



black hole masses in early-type galaxies at ^ 10** Mq (e.g., 
Shethetal. 2003). Total (integrated) quasar lifetimes esti- 
mated from observations are inferred to increase with increas- 
ing black hole mass as we predict (Yu & Tremaine 2002), and 
furthermore, the Eddington ratios of observed quasar samples 
are seen to increase systematically with redshift, as the sample 
becomes increasingly dominate d by luminosities above the 
break in the luminosity function jMcLure & Dunlot^2004^ ■ 

Moreover, observations show that the evolution of the lu- 
minosity function with decreasing redshift is driven by a de- 
crease in the characteristic mass scale of actively accreting 
black holes (e.g., |Heckman et al. 2004), which can be ex- 
plained in our model by the relation of the observed lumi- 
nosity function to the peak in the distribution of active black 
hole masses n(Lpeak)- This observation, however, has caused 
considerable confusion, as observations of both r adio-q uiet 
fe/oo & Urrv 2002) and radio-loud (O'Dow d etalJl2002l) lo- 
cal (low redshift) AGN indicate that nuclear and host lu- 
minosities are uncorrelated, implying that n uclear luminos- 
ity does not depend on black hole mass (Heckman et alj 
J004), and therefore that the primary variable determining 
the nuclear luminosity is the Eddington ratio, with the lumi- 
nosity func tion spanning a broad range in Eddington ratios 
SHao et al ]|2005). Furthermore, observations show that this 
is not true of high redshift quasars, as both direct estimates 
of acc retion rates (e.g., Vestergaard 2004; McLure & Dunlo^ 
120041) and the fact that their high luminosities would yield un- 
reasonably large black hole masses rule out substantially sub- 
Eddington accretion rates for most objects. Many previous 
empirical and semi-analytical models could not simultane- 
ously account for these observations. To explain just the low- 
redshift observations, such models adopt tunable distributions 
of Eddington ratios fitted to the data. However, both these 
observations are consequences of our interpretation of the lu- 
minosity function, as observations of local AGN and the low- 
redshift luminosity function are dominated by quasars below 
the break in the luminosity function, which are undergoing 
sub-Eddington growth and span a wide range of Eddington 
ratios, while observations at high redshift are dominated by 
bright objects at or above the break in the luminosity function, 
which are undergoing Eddington-limited (or near Eddington- 
limited) growth near their peak luminosity (see §|5j- 

(2) The second paradigm shift indicated by our modeling is 
that quasar obscuration is not a static or quasi-static geomet- 
ric effect, but is primarily an evolutionary effect. The physical 
reasoning for this is simple: the massive gas inflows required 
to fuel quasar activity produce large obscuring columns, and 
so column densities are correlated with quasar luminosity. 
The basic picture of buried quasar activity associated with the 
early growth of supermassive black holes and starburst activ- 
ity h as been proposed pr eviousl y arul studie d for some time 
(e.g.,'Sanders & Mirabel" 1996'; Fabian 1999), but our model- 
ing allows us to describe the evolution of obscuration in a self- 
consistent manner, defining obscured and unobscured phases 
appropriately and identifying dynamical correlations between 
the column density distribution and instantaneous and peak 
luminosities. 

There is substantial observational support for this picture. 
Point-like X-ray sources have been observed in many bright 
sub-millimeter or infrared and starburst sources, with es- 
sentially all very luminous infrared ga laxies showing ev- 
idence of buried quasar activity (e.g., Sanders & Mirabejl 
i996; Komo ssaeta l. 2003; Ptak et aL 2003 )- indicating si- 
multaneous buried black hole growth and star formation 
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at redshifts correspondin g to peak quasar activity (z > 1) 
jAlexander et alJl2005allj) . The buried black holes in high- 
z starbursting galaxies appear to be active but undermassive 
compa red to the quie scent galaxy black hole-stellar mass rela- 
tion (Borvs et al. 2005), implying that they are rapidly grow- 
ing in the starburst but have not yet reached their final masses, 
presumably set in the subsequent "blowout" phase. Simi- 
larly, observations suggest that obscured AGN are signifi- 
cantly more likely to exhibit strong sub-millimeter emission 
characteristic of star formation, implying both that obscured 
black hole growth and star formation are correlated and that 
obscuration mechanisms (responsible for re-radiation in the 
submm and IR) may be primarily isotropic in at least some 
cases (e.g.. Page et al. 2004; Stevens et al. 2005). Evidence 
from quasar emission line structure (e.g., Kuraszkiewi cz et al] 
12000: Tran 2003), directly related to the inner broad-line re- 
gion, suggests that isotropic obscuration of quasars can be 
important, in contradiction to angle-dependent mode ls. Fi- 
nally, many obse rvations (e.g. , Steffen et al. 2003; Ued a et al] 
2p03; Hasinger 2004; Grimes. RawHnes. & Willott 2004 
Sazon ov & Revnivtsev 2004; Bargeretal. 2005; Simpsog 
2005h indicate that the fraction of broad-line or obscured 
quasars is a function of luminosity, which cannot be 
accounted for in traditional static "torus" models (e.g., 
|Antonucci'1993) or reproduced even by modified luminosity- 
dependent torus models ( Lawrence 1 99 1 ), an observation that 
is explained by our model (see §|4]for a detailed discussion). 

Much of the obscuration in our modeling comes from large 
scales, arising from the inner regions of the host galaxy on 
scales ~ 50pc or larger While our resolution limits prevent 
our ruling out the possibility of gas collapse to a dense, ~pc 
scale torus surrounding the black hole, during the peak ob- 
scured phases of the final merger, our simulations indicate 
that these large scales dominate the contribution to the col- 
umn density, with quite large columns, which should be ob- 
servationally testable. Indeed, this is suggested by the typi- 
cal scales of obscuration in starbursting systems (e.g. Soifer 
et al. 1984a,b; Sanders et al. 1986, 1988a,b; for a review, 
see e.g. Soifer et al. 1987), given that, as discussed above, 
the dominant obscured phase of growth is closely associated 
with a starburst as implied observationally ( Alexander et a^ 
|l)05a b; B orvs et al.l2005l) . 

Observations of polarized light in intrinsically bright Type 
II AGN with unobscured luminosities typical of quasars (as 
opposed to local, dim Seyfert II objects in relaxed hosts) show 
scattering on large scales ~kpc, and in some cases obscura- 
tion clearly generated over scales extending beyond the host 
galaxy in the form of distortions, tidal tails, and streams from 
interactions and major mergers (Zakamska et al. 2004, 2005j). 
The angular structure seen in these observations is consistent 
with our mo deling . Moreover, in optically faint X-ray quasars 
(e.g. Donlev et al. 2005) it appears that obscuration is gener- 
ated by the host galaxies, and is directly related to host galaxy 
morphologies and line-of-sight distance through the host. The 
critical point is that, regardless of the angular structure of 
obscuration, typical column densities are strongly evolving 
functions of time, luminosity, and host system properties, and 
the observed distribution of column densities is dominated 
by these effects, not by differences in viewing angle across 
a uniform population. This is the case in our modeling as the 
lognormal dispersion (across different lines of sight) in col- 
umn densities is a^^ ^ 0.4 for a given simulation at some in- 
stant, whereas typical column densities across simulations, as 
a function of instantaneous and peak luminosities, span sev- 



eral orders of magnitude from A^h ^ 10"^ - 10^^ cm ^. 

8.2. Specific Predictions of our Model 
Our predictions include: 

• Quasar Lifetimes: We find that for a particular source, 
the quasar lifetime depends sensitively on luminos- 
ity, with the observed lifetime in addition depending 
on the observed waveband. Intrinsic quasar lifetimes 
vary from Iq ~ 10*'- 10*^ yrs, with observable lifetimes 
~ 10^ yrs in optical B-band cHopkins et. al... 2005aLh) . in 
good agr eement with o bservational estimates (for a re- 
view, see lMartinl2004 . 

• Luminosity Functions: Using a parameterization of the 
intrinsic distribution of peak luminosities (final quasar 
black hole masses) at a given redshift, our model of 
quasar lifetimes allows us to reproduce the observed 
luminosity function at all luminosities and redshifts 
z = 0-6. Although this is an empirical determination 
of the peak luminosity distribution, it im plies a new in- 
terpretation of the luminosity function (Hopkins et alj 
|2005c), which provides a physical basis for the ob- 
served "break" corresponding to the peak in the peak 
luminosity distribution. Moreover, the faint end slope 
is not determined by our empirical fitting procedure, 
but instead by the dependence of the quasar lifetime 
on luminosity, with its value and redshift evolution pre- 
dicted by our modeling ( Hopkins et al. 20051). The 
evolution of typical column densities in different stages 
of merger activity produces a significant population 
of obscured quasars, accounting for the difference be- 
tween hard X-ray (e.g., Ueda et al. 2003), soft X-ray 
(e.g., Mivaii et al. 2001), and optical B -ban d (e.g., 
ICroom et al„2004) luminosity functions (§ 13. 2t . 

• Column Density Distributions: The evolution of the 
column densities in our simulations reproduces the 
observed distribution of columns in optically-selected 
quasar samples, when the appropriate selection criteria 
are applied (Hopkins et al. 2005b), as well as complete 
column distributions in hard X-ray selected samples 
(§ 13. 3t . Column density evolution over the course of a 
merger yields a wider observed distribution of columns 
than that produced across different viewing angles at a 
given point in a merger 

• Broad Line Luminosity Function and Fraction: Using 
our simulations to estimate when quasars will be ob- 
servable as broad-line objects (either based on the ra- 
tio of quasar to host galaxy optical B-band luminos- 
ity or the obscuring column density), we reproduce 
the luminosity function of broad-line quasars in hard 
X-ray selected samples as well as optical broad-line 
quasar surveys, and the fraction of broad-line quasars 
in a given sample as a function of luminosity, to bet- 
ter precision than traditional or luminosity-dependent 
(but non-dynamical) torus models which are fitted to 
the data (§ 14. 2t . By providing an a priori prediction 
of the broad-line fraction as a function of luminosity 
and redshift which depends systematically on the typ- 
ical quasar host galaxy gas fraction, we propose that 
observations of the broad line fraction at different red- 
shifts can be used to constrain the gas fraction of quasar 
hosts and its evolution with redshift. 
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Active Black Hole Mass Functions: Using our pre- 
scription for deciding when objects will be visible as 
"broad-line" quasars, we predict the distribution of low- 
redshift, broad-line and non-broad line active quasar 
masses, in good agreement with observations from the 
SDSS, with expected incompleteness in the observed 
sample at low Mbh < 10^ M© black hole masses ('§ I4.3> . 
This is a new prediction which can be tested in greater 
detail by future observations, and our calculations allow 
us to model the differences in active black hole mass 
functions of the Type 1 and Type 11 populations. The 
width of the expected broad-line black hole mass func- 
tion depends significantly on the model of quasar life- 
times, enabling such measurements to probe the statis- 
tics of quasar evolution. 

Eddington Ratios: We determine Eddington ratio distri- 
butions from our simulations, given the peak luminosity 
distribution implied by the observed quasar luminosity 
function. The predicted distribution, once the appropri- 
ate observed magnitude limit is imposed, agrees well 
with observations at both low (z < 0.5) and high (1 .5 < 
z < 3.5) redshifts (§ |5Jl. As noted above, our inter- 
pretation of the luminosity function explains seemingly 
contradictory observations of Eddington ratios at dif- 
ferent redshifts. There is even a suggestion (Cao & x3 
|2005) that the evolution of quasars seen in our simula- 
tions (with bright phases in mergers and extended re- 
laxation after) can account for observa tions of bimodal 
Eddington ratio distributions at z ( Marchesini et alJ 
|2004), when coupled with an appropriate description 
of radiatively inefficient accretion phases, although it is 
possible that many of these low-redshift black holes are 
not fueled by mergers, especially in e.g. low-luminosity 
Seyferts. 

Relic Black Hole Mass Function: With our model for 
quasar lifetimes, the luminosity function at a given red- 
shift implies a birthrate of sources with given peak lu- 
minosities, n(Lpeak), which translates to a distribution 
in final black hole masses. Integrating this over red- 
shift, we predict the present-day mass distribution and 
total mass density of supermassive black holes. They 
agree well with observational estimates inferred from 
local populations of galaxy spheroids. In our picture, 
these spheroids are produced simultaneously with the 
supermassive black holes they harbor (§|6|l. We demon- 
strate that the integrated supermassive black hole den- 
sity, quasar flux density, and number counts in differ- 
ent wavebands can be reconciled with a radiative effi- 
ciency = 0.1, satisfying the constrain ts of counting 
arguments such as th at of i SoltanI il982ft . Further, we 
show in § 12.41 and § 14.11 that the corrections to such 
observational arguments based on optical quasar sam- 
ples are small (order unity) when we account for the 
luminosity dependence of quasar lifetimes, despite an 
extended obscured phase of quasar growth. In other 
words, although a quasar spends more time obscured 
than it does as a bright optical source, the total mass 
growth and radiated energy are dominated by the final 
"blowout" stage visible as a bright optical quasar. 

X-ray Background: The integrated quasar spectrum 
from our models of quasar lifetimes and column den- 
sities as a function of instantaneous and peak luminosi- 



ties can be combined with the birthrate of quasars with 
a given peak luminosity to give the integrated cosmic 
background in any frequency range. We predict both 
the normalization and shape of the X-ray background 
from ^ 1 - lOOkeV, with our modeling accounting for 
quasar obscuration as an evolutionary process (with a 
corresponding population of Compton-thick objects), 
avoiding any need for arbitrary assumptions about ad- 
ditional obscured populations (§ 17. 2> . For any model 
in which the quasar spectrum depends on luminosity or 
accretion rate, we demonstrate that a proper modeling 
of the quasar lifetime is critical to reproducing observed 
backgrounds. 

• Correlation Functions: In lLidz et alJ i2005l) . we predict 
the quasar correlation function and bias as a function 
of redshift and luminosity using our model, and com- 
pare it to that expected using "light bulb" or exponential 
light curves. As most quasars in our modeling have a 
characteristic peak luminosity or final black hole mass 
corresponding to the peak of the n(Lpeak) distribution, 
they reside in hosts of similar mass, and there is lit- 
tle change in bias with luminosity at a given redshift, 
in contrast to idealized models for the quasar lifetime 
and luminosity function. Our predicted bias agrees well 
with the observations of Croom et al. (2005), who also 
find no evidence for a dependence of the correlation 
on quasar luminosity at a given re dshift, as we ex pect. 
In fact. IPor ciani. Magliocchetti, & Norberd J2004I) and 
ICroom etal.1 (i2005i,) find that their observations can be 
explained if quasars lie in hosts with a constant char- 
acteristic mass ^ 2 X lO'^M© {h = 0.7). If we consider 
their redshift range 1 - 2, we predict the quasar pop- 
ulation will be dominated by sources with Lpeak = i*(z), 

which given MR„(LDeak) and using theMBH-Mhaio rela- 
tion of Wvithe & Loeb (2003) yields a nearly constant 
characteristic host hal o mass ^ 1 -2 x lO'^M©, in good 
agreement. Similarlv. lAdelberger & Steidell ( 12005 ) find 
that the quasar-galaxy cross-correlation function does 
not vary with luminosity, implying with ^ 90% confi- 
dence that faint and bright quasars reside in halos with 
similar masses and that fainter AGN are longer lived, 
strongly disfavoring traditional "light bulb" and expo- 
nential light curve models. Furthermore, Hennawi et aQ 
([2005) find an order of magnitude excess in quasar clus- 
tering at small scales < 40/1"' kpc, with the correlation 
function becoming progressively steeper at sub-Mpc 
scales, suggesting that quasar activity is triggered by 
interactions and mergers. 

• Host Galaxy Properties: Because black hole growth 
and spheroid formation occur together in our picture, 
our modeling allows us to describe relationships be- 
tween black hole and galaxy properties. For exam- 

Sile, we reproduce both the observed Mbh - cr relation 
Di Matteo et al. 2005; Robertson et al. 2005b) and the 
fundamental plane of elliptical galaxies (Robertson et 
al., in preparation). Since we also reproduce the distri- 
bution of relic black holes inferred from the z = Q dis- 
tribution of spheroid velocity dispersions or luminosity 
functions using the observed versions of these relations, 
our match to these relations indicates that we also re- 
produce these distributions of host spheroid prope rties. 
We consider this in detail in lHopkinsetaf](l2005el) . and 
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find that we are able to account for a wide range of host 
galaxy properties, including luminosity and mass func- 
tions, color-magnitude relations, mass-to-light ratios, 
and ages as a function of size, mass, and redshift. With 
our modeling of the quasar lifetime as motivated by our 
simulations, the evolution and distribution of properties 
of red-sequence galaxies and the quasar population are 
shown to be self-consistent, which is not the case for 
idealized models of quasar evolution. 

Aside from an empirical estimate of the distribution of peak 
quasar luminosities «(Lpeak)> we determine all of the quantities 
summarized above self-consistently from the input physics of 
our simulations, including a physically motivated dynamic ac- 
cretion and feedback model in which black holes accrete at the 
Bondi rate determined from the surrounding gas, and 5% of 
the radiant energy couples thermally to that gas. Beyond this, 
our simulations enable us to calculate the various predictions 
above a priori, without the need for additional assumptions or 
tunable parameters. 

We compare each of these predictions to those obtained us- 
ing idealized descriptions of the quasar lifetime, i.e. "light- 
bulb" and exponential light curve (constant Eddington ratio) 
models, and the column density distribution, i.e. standard and 
"receding" (luminosity-dependent) torus models. We fit all 
these (along with our full model) to the observed luminos- 
ity function in the same manner (allowing the same degree 
of freedom to ensure that they all yield the same observed 
luminosity function), and we fit the free parameters of these 
tunable models (e.g. typical Eddington ratios and quasar life- 
times for the "light-bulb" or exponential models, typical col- 
umn densities and torus scalings for the torus models) inde- 
pendently to each observation to maximize their ability to re- 
produce observations. However, we still find better agreement 
between our model (with no parameters tuned to match obser- 
vations) and the observations in nearly every case where the 
tunable phenomenological model is not guaranteed to repro- 
duce the observation by construction. The one exception is 
the relic supermassive black hole mass function, for which 
the predictions of our modeling and idealized lifetime models 
are essentially identical, reflecting the fact that in both cases 
black hole growth is dominated by bright, optically observ- 
able, high Eddington ratio phases. 

Moreover, the best-fit parameters for the idealized mod- 
els, when fitted independently to each observation, are not 
self-consistent. For example, calculations of the black hole 
mass function imply high Eddington ratios / ^ 0.5 - 1 (e.g., 
^u & Tremaine 2002), and our fit to the active black hole 
mass function ( Hec kman et al. 2004) suggests / ^ 1, but the 
observed distribution of accretion shows a typical / ^ 0.3 
(IVestergaard 

20(3, and fitting to the broad-line fraction as 
a function of luminosity with our full obscuration model but 
these lifetime models implies a lower / ^ 0.05. Likewise, 
fitting torus models to the X-ray background suggests typ- 
ical column densities through the torus of A^h ~ 10^^ cm"^ 
(e.g.,|Treister & Urrv 2005), while fitting to the observed col- 
umn density distributions (.Treister et al.l2004l iMainieri et al.l 
1^05) suggests equatorial columns A^h ^ lO^'^cm"^. Clearly 
then, reproducing the observations listed above, and in partic- 
ular doing so self-consistently, is not implicit in any model 
which successfully reproduces the quasar luminosity func- 
tion, even at multiple frequencies. 

8.3. Further Testable Predictions of our Model 



Our model for quasar evolution makes a number of obser- 
vationally testable predictions: 

• Quasar lifeti mes are only w eakly constrained by obser- 
vations (e.g. ' Martiniir2004ft . but future studies may be 
able to measure both the hfetime of individual quasars 
and the statistical lifetimes of quasar populations as 
a function of luminosity. We describe in detail our 
predictions for the evolution of individual quasars and 
quantify their lifetimes in § |2l and further predict the 
distribution of both integrated and differential lifetimes 
in an observed sample as a function of luminosity. This 
should provide a basis for comparison with a wide 
range of observations, with the most important predic- 
tion being that the quasar lifetime should increase with 
decreasing luminosity. 

• For a reasonably complete, optically selected sample 
above some luminosity, the distribution of observed 
column densities should broaden to both larger and 
smaller A^h values as the minimum observed luminosity 
is decreased, as both intrinsically faint periods with low 
column density and intrinsically bright periods with 
high column density become observable. 

• Similarly, the Eddington ratio distribution should be a 
function of observed luminosity, with a broad distribu- 
tion of Eddington ratios down to / ~ 0.01 -0.1 at lumi- 
nosities well below the break in the observed luminos- 
ity function, and a more strongly peaked distribution 
about I ^ 0.2- 1 for luminosities above the break (Fig- 
urel20l. 

• In our interpretation, the bright and faint ends of the 
luminosity function correspond statistically to similar 
mixes of galaxies, but in various stages of evolution; 
whereas in all other competing scenarios, the quasar 
luminosity is directly related to the mass of the host 
galaxy. Therefore, any observational probe that differ- 
entiates quasars based on their host galaxy properties 
such as, for example, the dependence of the cluster- 
ing of quasars on luminosity, or the host stellar mass 
and size as a function of luminosity (although we cau- 
tion that this is somewhat dependent of the modeling 
of star formation in mergers), can be used to discrim- 
inate our picture from older models. We present a 
detailed prediction of the quasar correlation function 
based on our modeling for comparison with observa- 
tions in Lidz et al. (2005). 

• Our distribution «(Lpeak) directly translates to a black 
hole merger rate, as a function of mass, in our model- 
ing, allowing a detailed prediction of the gravitational 
wave signal from black hole mergers as a function of 
redshift. 

• The broad line fraction as a function of luminosity, de- 
fined by requiring that "broad-line" objects have an ob- 
served B-band luminosity above a fraction /bl of that 
of their host galaxy, is a prediction of our model quasar 
and galaxy light curves. However, the uncertainties are 
large, primarily because different observational sam- 
ples have varying sensitivity to quasar vs. host galaxy 
optical light. Furthermore, the host galaxy gas fraction 
and /bl are degenerate in these predictions - a well- 
defined observational sample complete to some /bl can 
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constrain our modeling of quasar fueling and the rela- 
tion between quasar and host galaxy light curves. In 
particular, such observations, either by measuring the 
faint-end shape of the "broad-line" quasar luminosity 
function or the mean "broad-line" fraction at a given lu- 
minosity as a function of redshift, can constrain the gas 
fractions of quasar host galaxies and their evolution, es- 
sentially a free parameter in our empirical modeling. 

• We also predict the distribution of active, low-redshift 
black hole masses in §0] These predictions can be com- 
pared to mass functions for active black holes from nu- 
merous quasar surveys, which should include improved 
mass functions of the entire quasar population complete 
to lower luminosities as well as future mass functions 
for the population of active broad-line AGN. We pro- 
vide predictions for the black hole mass function of all 
active quasars, and for just the "broad-line" population 
(as a function of the survey selection). 

• Because the evolution of spheroids and supermassive 
black holes is linked in our modeling, with each affect- 
ing the evolution of the other, we can also use the distri- 
bution of observed quasar properties to predict galaxy 
properties such as number counts, spheroid masses and 
luminosities, and colors as a function of redshift. For 
the calculation and discussion of these predictions, see 
[Hopkins et al. (2005e). 

• In our model, the growth of supermassive black holes 
is dominated by galaxy mergers. Therefore, at any 
given redshift, the mass (and as a consequence, lu- 
minosity) function of galaxy mergers should have a 
similar shape to our distribution of quasar birthrates, 
n(Lpeak), distinct from the shapes of either the quasar 
or total galaxy luminosity functions. Indeed, prelim- 
inary observational estimates of both the me rger lu- 
mino s ity function (e.g., Xu et al. 2004; Conseli ce et alJ 
120031 iWolf et alJ 12005) and quasar host gala xy lu- 
mino sity function iBahcall et al. 1997; Hamiltor i et alJ 
12002'), primarily at low redshifts, appear be consistent 
with this expectation. Theoretically, it may be possi- 
ble to predict the merger luminosity function using ei- 
ther cosmological simulations or semi-analytical mod- 
els; we discuss this further in §|5] 

8.4. Mock Quasar Catalogs 

In principle, our modeling can be used to predict the distri- 
butions of quasar luminosities in various wavebands, column 
densities, active black hole masses, and peak luminosities for 
a wide range of observational samples, but it is impractical 
for us to plot predictions of these quantities for all possible 
sample selection criteria. To enable comparison with a wider 
range of observations, we have used our modeling and the 
conditional probability distributions for these quantities from 
our simulations to generate Monte Carlo realizations of quasar 
populations, which we provide publicly via ftp'*. 

At a particular redshift, we use our fitted n(Lpeak) distribu- 
tion and our suite of simulations to generate a random popula- 
tion of mock "quasars." We first generate the peak luminosi- 
ties of each "quasar" according to the fitted n(Lpeak) at that 
redshift. For each object, we then use the probability of be- 
ing at a given instantaneous luminosity in simulations with a 

'* |f tp : //cf a-f tp ■ harvard. edu/pub/phopkins/qso_catalogs /I 



similar peak luminosity to generate a current bolometric lumi- 
nosity. In practice, we calculate the P{L |Lpeak) distribution by 
summing w(Lpeak, ipeak, i) X P(L |Lpeak. i), where Lpeak is the 
mock quasar peak luminosity, Lpeak, i is the peak luminosity 
of each simulation and w(Lpeak, Lpeak, i) is a Gaussian weight- 
ing factor (oc exp(-log^(Lpeak/Lpeak, i)/2(0.05)2)). Knowing 
the instantaneous bolometric luminosity L and peak luminos- 
ity Lpeak, we then follow an identical procedure to determine 
the joint distribution P{X \ L, Lpeak) of each subsequent quan- 
tity X, from simulations with similar L and Lpeak- We have 
compared this with Monte Carlo realizations based on our fit- 
ted probability distributions in this paper, and find that essen- 
tially identical results are achieved for e.g. the distribution of 
L and Lpeak, and column densities in phases of growth not near 
peak luminosity. However, this modeling is not identical for 
e.g. the distribution of Eddington ratios and column densities 
around L ~ Lpeak, which reflects the fact that our fits to the 
Eddington ratio distribution (§ |5) are rough and that our fits 
to the column density distribution do not apply to the final 
"blowout" phase of quasar evolution (as discussed in detail in 

§13. 

For each mock quasar, we generate a peak luminosity, fi- 
nal (post-merger) black hole mass, instantaneous bolomet- 
ric luminosity, intrinsic (un-attenuated) B-band (uLi, ?X v = 
4400A), soft X-ray (0.5-2 keV), and hard X-ray (2-10 keV) 
luminosity, observed (attenuated using the generated column 
density and the reddening/dust extinction modeling described 
in § 12.21 with SMC-like reddening curves and extinction fol- 
lowing e.g. Pei 1992, Morrison & McCammon 1983) B-band, 
soft X-ray, and hard X-ray luminosities, column density of 
neutral hydrogen, column density of neutral-nionized hydro- 
gen, and instantaneous black hole mass. The intrinsic lumi- 
nosities in each band are calculated using the bolometric cor- 
rections described in Marconi et al. ( 2004), which account for 
the luminosity dependence of t he op tical-to-X-ray luminos- 
ity ratio aox (as discussed in § I3.2l i. and then attenuated to 
give the observed luminosities. We also provide intrinsic and 
attenuated luminosities in each waveband using the constant 
bolometric corrections of Elvis et al. ( 1994), but we caution 
that these are not calculated in a completely self-consistent 
manner, as our assumed bolometric luminosity function to 
which we fit the «(Lpeak) distribution is based on using the 
luminosity-dependent bolometric corrections. We do not di- 
rectly calculate Eddington ratios, as these are defined differ- 
ently in many observed samples, but they should be calculable 
with the given luminosities and black hole masses. 

We calculate these quantities for a mock sample of ~ lO' 
quasars at each redshift z = 0.2, 0.5, 1, 2, and 3. Most of 
these quasars are at luminosities orders of magnitude below 
those observed, therefore for space considerations and be- 
cause our predictions become uncertain at low luminosities, 
we retain only the 10^ quasars with brightest bolometric lu- 
minosities at each redshift. This introduces some uncertainty 
in our statistics at the lowest luminosities in any given band, 
but these luminosities are generally still well below those ob- 
served in most samples. At any luminosity, but especially at 
the brightest luminosities, there is also a significant amount 
of effective "noise" owing to our incomplete sampling of the 
enormous parameter space of possible mergers, and decreas- 
ing total time across simulations spent at large luminosities, 
which can be estimated from e.g. Figures |8] and Finally, 
at each redshift, we generate two distributions, reflecting the 
~ Icr range in «(Lpeak), and roughly parameterizing the de- 
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generacies in our fit to the observed luminosity functions and 
uncertainty in the faint-end of n(Lpeak) - "Fit 1" has a lower 
(lower peak in «(Lpeak)), with a larger a, (broader n'(Lpeak) 
distribution), and "Fit 2" has a higher L» and smaller cr» (more 
narrowly peaked n(Lpeak) distribution). We show a few exam- 
ple "quasars" from our z = 0.2 mock distribution in Tabled 
to demonstrate the format and units used. 

8.5. Starburst Galaxies 

Although we do not yet model the re-radiation of absorbed 
light by dust or the contribution of stellar light to quasar host 
IR luminosities, including these in our picture for quasar evo- 
lution will enable us to predict luminosity functions in the 
IR and sub-mm and their evolution with redshift. We can at 
this point, however, estimate if our model for quasar lifetimes 
and merger-driven evolution with «(Lpeak) is consistent with 
the observed distribution of ultraluminous infrared galaxies. 
Naively, we might expect that since the obscured quasar phase 
has a duration up to ~ 10 times that of the optically observable 
quasar phase, there should be ~ 10 times as many ULIRGs as 
bright optical QSOs. But, this neglects the complicated, lumi- 
nosity dependent nature of quasar lifetimes. 

Given that the bright quasars we simulate attain, during 
their peak growth phase, an intrinsic luminosity comparable 
to that of the host starburst, and that this period of peak growth 
has a similar duration to the starburst phase (se e Figure [T3l 
and'Di Matteo et al.""2005; Sprinael et al."2005b'), we can es- 
timate (roughly) the ULIRG bolometric luminosity function 
from our bolometric quasar luminosity function. Thus, the 
more accurate comparison to the ULIRG luminosity func- 
tion is with the hard X-ray quasar luminosity function, as 
this recovers (and at some luminosities can be dominated by) 
"buried" quasars in starburst phases. This is only applica- 
ble above the break in the luminosity function, where quasars 
are undergoing peak quasar growth. Below the break, quasars 
are, on average, sub-Eddington and can have luminosities well 
below that of their star-forming hosts (see Figure [O}, so we 
expect our quasar luminosity function to be significantly shal- 
lower than the ULIRG luminosity function at these luminosi- 
ties. Note also that this does not imply that ULIRGs are all 
AGN-dominated, as the starburst and peak AGN activity can 
be (and generally are) somewhat offset, but only says that the 
lifetime curves at the bright end should be similar 

Considering the luminosity function at z = 0.15, then, 
we expect ULIRG densities ii$/dMboi ^ 3 x 10"^ and 9 x 
10"** Mpc-3 mag-i at L - 1 .6 x IO'^Lq and 2.5 x IO'^Lq, re- 
spectively. These estimates are consistent with the observed 
density in the IRAS 1 Jy Survey (Kim & Sanders, .1998:.) at 
a mean redshift z = 0.15, with d$/dMboi ~ 5 x 10"^, 7 x 
10"** Mpc'-'mag"' (rescaled to our cosmology), and as ex- 
pected, our quasar luminosity function slope becomes signif- 
icantly shallower than the observed 1 Jy survey luminosity 
function slope below L ^ lO^'-lO'^L©, roughly the break 
luminosity of the luminosity function. We predict these den- 
sities to change with redshift according to the evolution of 
n(Lpeak), decreasing by a factor ^ 1.5 at z = 0.04, in good 
agreement with the evolution of IRAS ULIRG luminosity 
functions ( Kim & Sanderslfl998h . Likewise, at z ~ 1 - 3, 
we predict a mean space density <I>(L > lO'^L©) ~ 1-3 x 
lO'' Mpc"^, in agreement with the ~ 5 x 10^ Mpc"^ density of 
such sources expected to reproduce the observed cumulative 
source density 4x1 0^ deg~^ of 1 mJy 850 /im SCUBA sources 
dBargeret al.lll999h . Furthermore, our prediction of the frac- 
tion of buried AGN and its evolution with redshift agrees well 



with determinations from X-ray samples dBarger et al.ll200^ 
and recent Spitzer results in the mid and near-infrared at z ~ 2 
dMartinez-Sansigre et alJ2005h . 

8.6. AGN not Triggered by Mergers 

Some low redshift quasars (e.g. Bahcall et al. 1996) and 
many nearby, low-luminosity Seyferts appear to reside in or- 
dinary, relatively undisturbed galaxies. Our picture for quasar 
evolution does not immediately account for these objects be- 
cause we suppose that nuclear activity is mainly triggered by 
tidal torques during a merger 

This work is primarily concerned with the origin of the ma- 
jority of the mass in spheroids and supermassive black holes, 
and as a consequence, the relation of this to the abundance and 
evolution of quasars and the cosmic X-ray background. Based 
on our present analysis, we believe that a merger-driven pic- 
ture can account for the main part of each of these, and, as 
described earlier, that the most relevant phase in the history 
of the Universe to these phenomena appears to have been at 
moderate redshifts, z ^ 2.5 to z ^ 0.5. 

Our model does not exclude other mechanisms for trigger- 
ing AGN and it is likely that a variety of stochastic or con- 
tinuous processes are relevant to nuclear activity in undis- 
turbed disks and residual low-level accretion in relaxed sys- 
tems. This is not contrary to our picture because most of the 
total black hole mass density in the Universe is in spheroid- 
dominated systems. The principal requirement of our model 
is that AGN activity in undisturbed galaxies should not con- 
tribute a large fraction of the black hole mass density in 
the Universe, to avoid spoiling tight correlations between the 
black hole and host galaxy properties and producing too large 
a present-day black hole mass density in violation of the 
Soltan (1982) constraint. 

For example, if a molecular cloud passed through the cen- 
ter of our Galaxy near Sgr A*, it is possible that the Milky 
Way would resemble a Seyfert for some period of time. 
Alternatively, it has long been recognized that mass loss 
from normal stellar evolution of bulge stars or stellar clus- 
ters near the centers of galaxies can provi de a continuous 
supply of fuel for low-level accretion (e.g.. iMcMillan et al.l 
1981; MacDonald & Bailev 1981; Shull 1983). Typical galac- 
tic stellar mass loss rates (M - IMoyr"' (IO^Mq)"') yield 
Bondi-Hoyle accretion rates ~ 1 0~^ - 1 0"'* of Eddington in re- 
laxed, dynamically hot systems; and mass loss rates from O 
and W-R stars (M ~ \0~^Mq yr"' (lOM©)"') in young, dense 
star clusters near the centers of galaxies with sufficient cold 
gas for continued star formation can yield rates as high as 
~ 10"^ of Eddington. 

Even though these fueling mechanisms do not involve 
mergers, the scenario we have discussed might still be rele- 
vant to the origin of these black holes. Of course, the black 
holes and spheroids in disk-dominated systems may have pro- 
duced in a manner that did not involve mergers. Alternatively, 
most of the black hole mass in these objects (which is small 
compared to that in spheroid-dominated galaxies) could have 
been assembled long ago in mergers with bright quasar phases 
and then these "dead" quasars are resurrected sporadically by 
other fueling mechanisms. 

Independent of how these black holes were formed, ele- 
ments of our modeling may still account for certain observed 
properties of Seyferts. The observed Seyfert luminosity func- 
tion appears to join smoothly onto the quasar luminosity func- 
tion (,Hao et al. 2005 ). It is not obvious that this would be 
the case if the two types of objects are produced by entirely 
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distinct mechanisms. In addressing this, it is useful to sepa- 
rate the process by which gas is delivered to the black hole 
from the subsequent evolution that determines the observed 
activity. In our picture, gas is delivered to the black hole by 
gravitational torques during a merger, but other mechanisms, 
like bar-induced fueling may be important for objects such as 
Seyferts. Regardless, the induced activity may be generic, if 
black hole growth is self-regulated in the way we describe it 
in our simulations. 

In Hopkins et al. (2005f) we show that the faint end slope 
of the quasar luminosity function in our model is partly de- 
termined by the time dependence of the "blowout" phase of 
black hole growth. We derive an analytical model for this us- 
ing a Sedov-Taylor type analysis and show that the impact of 
this feedback depends on the mass of the host. This analy- 
sis does not depend on the fueling mechanism, only on the 
subsequent evolution. If this self-regulated growth applies to 
Seyferts as well (for example if Seyfert growth is regulated 
by a balance between accretion feedback and the spheroid po- 
tential, as expected if these objects obey a similar Mbh - o" 
relation), we would expect the Seyfert luminosity function to 
smoothly join onto the quasar one, even if the fuel is delivered 
in a different manner. 

9. CONCLUSIONS 

We have studied the evolution of quasars in simulations of 
galaxy mergers spanning a wi de region of parameter space. In 
agreement with earlier work ( Hopkins et al.ll2005a l). we find 
that the lifetime of a particular source depends on luminosity 
and increases at lower luminosities, and that quasar obscu- 
ration is time-dependent. Our new, large set of simulations 
shows that the lifetime and obscuration can be expressed in 
terms of the instantaneous and peak luminosities of a quasar 
and that these descriptions are robust, with no systematic de- 
pendence on simulation parameters. We have combined these 
results with a semi-empirical method to describe the cosmo- 
logical distribution of quasar properties, allowing us to predict 
a large number of observables as a function of e.g. luminosity 
and redshift. This approach also makes it possible to compare 
our picture to simpler models for quasar lifetimes and obscu- 
ration. 

In the model we examine, quasars are triggered by mergers 
of gas-rich galaxies, which produce inflows of gas through 
gravitational torquing, fueling starbursts and rapid black hole 
growth. The large gas densities obscure the central black hole 
at optical wavelengths until feedback energy from the growth 
of the black hole ejects gas and rapidly slows further accre- 
tion ("blowout"). Quasar lifetimes and light curves are non- 
trivial, with strong accretion activity during first passage of 
the merging galaxies and extended quiescent (sub-Eddington) 
phases leading into and out of the phase of peak quasar activ- 
ity associated with the final merger. The "blowout" phase in 
which the quasar is visible as a bright, near-Eddington optical 
source has a lifetime related to the dynamical time in the in- 
ner regions of the merging galaxies, which characterizes the 
timescale over which obscuring gas and dust are expelled, but 
the quasar spends a longer time at lower luminosities before 
and after this stage. These evolutionary processes have im- 
portant consequences which cannot be captured in models of 
pure exponential or "on/off" quasar growth. 

Our work emphasizes several goals for quasar and galaxy 
observations and theory. Observationally, it is important to 
constrain the faint end of the peak luminosity distribution; 
i.e. the low-mass active black hole distribution. Unfortu- 



nately, our modeling of quasar lifetimes implies that the faint- 
end quasar luminosity function is dominated by quasars with 
peak luminosities around the break in the luminosity func- 
tion, and can provide only weak constraints on the faint-end 
Lpeak distribution. However, there is still hope, as for ex- 
ample broad-line quasar activity is more closely associated 
with near-peak luminosities, and thus probing the faint-end 
of broad-line luminosity functions may in particular improve 
the estimates. Moreover, studies of the black hole mass dis- 
tribution (or the distribution of galaxy spheroids) as a func- 
tion of redshift, extending to small spheroid masses/velocity 
dispersions probes the faint end of n(Lpeak)- Other tech- 
niques, such as studies of faint radio sources at high red- 
shift ( Haiman, Ouataert, & Bower 2004) can similarly con- 
strain these populations. Furthermore, the calculations in this 
paper can be combined to better determine n(Lpeak), as, given 
a model for the quasar lifetime and obscuration, they all derive 
from this fundamental quantity. Additional observational tests 
of the modeling we have presented will provide an important 
means of constraining models for AGN accretion and feed- 
back; for example, the faint-end slope of the quasar lifetime 
depends on how the "blowout" phase occurs and could pro- 
vide a sensitive probe of feedback models, enabling the adop- 
tion of more realistic and sophisticated feedback prescriptions 
than we have thus far employed. Of course, improved con- 
straints on the luminosity function at all luminosities at high 
redshift remains a valuable means of testing theories of quasar 
evolution. 

Our simulations are based on isolated galaxy mergers, and 
thus do not provide a cosmological prediction for the distri- 
bution of peak luminosities n(Lpeak), merger rates, or mass 
functions - we instead have adopted a semi-empirical model, 
in which we use our modeling of quasar evolution to deter- 
mine these distributions from the observed luminosity func- 
tion. While this allows us to predict a large number of ob- 
servables and to demonstrate that a wide range of quasar and 
galaxy properties are self-consistent in a model of merger- 
driven quasar activity with realistic quasar lifetimes, future 
theoretical work in these areas should predict the distribution 
of peak luminosities n(Lpeak) and its evolution with redshift. 
These quantities are to be distinguished from the distribution 
of observed luminosities, as the two are not trivially related in 
our model or any other with a non-trivial quasar lifetime. 

Although the quasar birthrate as a function of peak lumi- 
nosity will be, in general, a complicated function of galaxy 
merger rates, gas fractions, morphologies, and other factors, 
we have parameterized it for comparison with the results of 
future cosmological simulations and semi-analytical models. 
This distribution is particularly valuable as a theoretical quan- 
tity because it is more directly related to physical galaxy prop- 
erties than even the complete (intrinsic) luminosity function, 
and additionally because theoretical modeling which success- 
fully reproduces this «(Lpeak) distribution is guaranteed to re- 
produce the large number of observable quantities we have 
discussed in detail in this work. We cannot determine the cos- 
mological context in our detailed simulations of the relatively 
small-scale physics of galaxy mergers, and conversely, cos- 
mological simulations and semi-analytical models cannot re- 
solve the detailed physics driving quasar activity in mergers. 
However, our determination of quasar evolution as a function 
of peak luminosity or final black hole mass can be grafted 
onto these cosmological models to greatly increase the effec- 
tive dynamic range of such calculations. Combined with our 
modeling, this would remove the one significant empirical el- 
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ement we have adopted, and allow for a complete prediction 
of the above quantities from a single theoretical framework. 

In these efforts, we emphasize that the mergers relevant to 
our picture are of a specific type. First, the merging galaxies 
must contain a supply of cold gas in a rotationally supported 
disk. Hot, diffuse gas, as in the halos of elliptical galaxies, 
will not be subject to the gravitational torques which drive 
gas into galaxy centers and fuel black hole growth. Clearly, 
gas-poor mergers are also not important for this process. Sec- 
ond, the mergers will likely involve galaxies of comparable, 
although not necessarily equal, mass, so that the gravitational 
torques excited are strong enough and penetrate deep enough 
into galaxy centers to drive substantial inflows of gas. The 
precise requirement for the mass ratio is somewhat ill-defined 
because it also depends on the orbit geometry, but mergers 
with a mass ratio larger than 10 : 1 are probably not generally 
important to our model. Simulations of minor mergers involv- 
ing galaxies with mass ratios < 10 : 1 (e.g. Hernquist 1989; 
Hernquist & Mihos 1995) have shown that for particular or- 
bital geometries, these events can produce starbursts similar 
to those in major mergers, leaving behind disturbed remnants 
with dynamically heated disks (e.g. Quinn et al. 1993; Mihos 
et al. 1995; Walker et al. 1996). It is of interest to estab- 
lish whether black hole growth can also be triggered in minor 
mergers, as these events may be relevant to weak AGN activ- 
ity like that in some Seyfert galaxies or LINERs. 

In summary, the work presented here supports the conjec- 
ture that many aspects of galaxy formation and evolution can 
be understood in terms of the "cosmic cycle" in Figure To 
be sure, much of what is summarized in Figure^has been pro- 
posed elsewhere, either in the context of observations or theo- 
retical models. Our modeling of galaxy formation and evolu- 
tion emphasizes the possibility that supermassive black holes 
could be responsible for much of what goes on in shaping 
galaxies, rather than being bystanders, closing the loop in Fig- 
ure [0 In this sense, black holes may be the "prime movers" 
driving galaxy evolution, as has been proposed earlier for ex- 



tragalactic radio sources (e.g. Begelman, Blandford & Rees 
1984; Rees 1984). It may seem counterintuitive that com- 
pact objects with masses much smaller than those of galaxies 
could have such an impact, but it is precisely the concentrated 
nature of matter in black holes that makes this idea plausible. 

Consider a black hole of mass Mbh at the center of a spher- 
ical galaxy of mass Mjph with a characteristic velocity disper- 
sion a. The energy available to affect the galaxy through the 
growth of the black hole will be some small fraction of its 
rest-mass, Zifeed ^ e/A^bhc^- This can be compared with the 
binding energy of the galaxy, Zibind ^ A^sphf^- Observations 
indicate that Mbh and Mjph are correlated and that, roughly 
Mbh (0.002 -0.005)M.sph (Magorrian et al. 1998; Mai'coni 
& Hunt 2003). Therefore, the ratio of the feedback energy to 
the binding energy of the galaxy is -Efeed/^sph > 10e/.-2O'3gQ, 
for an assumed efficiency of 1%, e/.-2 = e/0.01 and scaling 
the velocity dispersion to (J300 = cr/300 km/sec, as for rel- 
atively massive galaxies. This result demonstrates that the 
supermassive black holes in the centers of spheroidal galax- 
ies are by far the largest supply of potential energy in these 
objects, exceeding even the galaxy binding energy. When 
viewed in this way, if even a small fraction of the black hole 
radiant energy can couple to the surrounding ISM, then black 
hole growth is not an implausible mechanism for regulating 
galaxy formation and evolution; in fact, it appears almost in- 
evitable that it should play this role. 
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TABLE 1 

Mock Quasar Distribution Examples 



1 J 

^peak 




3l 


■*Mbh 








8,i 


^HX 


10 robs 


1 1 7 obs 
^SX 


12 robs 
^HX 


10.6 
10.4 


6.1 
6.4 


8.5 
8.7 


6.1 
6.0 


20.5 
22.2 


20.1 
22.0 


7.2 7.5 
7.4 7.6 


lA 6.8 
7.6 7.0 


7.6 7.0 
7.8 7.2 


7.2 7.4 
5.8 6.0 


7.4 6.8 

7.5 6.9 


7.6 7.0 
7.8 7.2 



'Peak quasar bolometric luminosity, logiQ(Lpeak/^'0) 
^Final (post-merger) black hole mass, logjo(Mgj^/M0) 

■^CuiTent (at time of "observation") intrinsic (no attenuation) bolometric luminosity, logjQ(L/LQ) 
'^Cun'ent black hole mass, logjQ(MBH/-'^0) 

^Total (neutral+ionized) hydrogen column density along the "observed" sighdine, logjQ(A'H/cm~^) 

^Neutral hydrogen column density along the "observed" sightline, logjQ(A^Hi/cm~") 

^Intrinsic (no attenuation) B-band luminosity, logjQ(Lg/LQ), where Lb = i^bLuj^ at = 4400A. 

Calculated with the luminosity-dependent bolometric con'ections from (Marconi et al. 2004; left), 

and constant (luminosity-independent) L = Il.SLe (Elvis et al. 1994; right). 
^Intrinsic soft X-ray (0.5-2 keV) luminosity, \og^Q(V^^ /Lq). 

Calculated with the luminosity -dependent bolometric corrections from (Marconi et al. 2004; left), 

and constant (luminosity-independent) L = 52.5Lsx (Elvis et al. 1994; right). 
^Intrinsic hard X-ray (2-10 keV) luminosity, \og^Q(L'fj^/LQ). 

Calculated with the luminosity-dependent bolometric corrections from (Marconi et al. 2004; left), 

and constant (luminosity-independent) L = 35.0Lhx (Elvis et al. 1994; right). 
^^"Observed" (with attenuation) B-band luminosity, logjQ(Lg''^/L0). 

Left and right use luminosity-dependent and luminosity-independent bolometric corrections, respectively, as L^. 
^^"Observed" soft X-ray luminosity, logjQ(L^^Y^0)- 

Left and right use luminosity-dependent and luminosity-independent bolometric conections, respectively, as . 

"Observed" hard X-ray luminosity, log]o(^^Jx /^o)- 

Left and right use luminosity-dependent and luminosity-independent bolometric conections, respectively, as L^ffx- 

The complete tables can be downloaded a t,f tp : / /cf a- f tp ■ harvard - ed'u/pub/phopkins/qso_catalogs/l 



